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Abstract 

We present results of numerical analysis of several simple models for the microstructure 



of a double auction market without intermediaries which were introduced in ||10|| . These 
markets can be represented as a set of buyers and a set of sellers, whose numbers vary 
in time, and who diffuse in price space and interact through an annihilation interaction. 
In this paper two models suggested in [jnj are studied - the minimal model and the two- 
liquid model. We perform computer simulations of the minimal model in order to verify 
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three of the liquidity scaling laws postulated in [|H|] . It is found that midmarket variance, 
bid-offer spread, and fluctuation of the bid-offer spread scale according to Dj J where D is 
the diffusion coefficient (trader volatility) and J is the deal rate. A logarithmic correction 
to the scaling law for midmarket variance is observed, but not for bid-offer spread or its 
fluctuation, because they are fundamentally different quantities. Scaling parameters are 
obtained. We show both analytically and numerically that the total number of traders in 
the market scales as JL 2 /4D, where L is the width of a price space. Time to midmarket 
sale (Ys) is found to scale as 1/J while its fluctuation goes as 0.73/ J. A "reduced" time 
(^reduced) is also studied, and found to scale in a non-trivial way. Asymmetric fluxes 
are introduced to the minimal model and analytical result derived in ||10|| , for the speed 
of the moving midmarket agrees with numerical results. Simulation of the two-liquid 
model which describes a market with both market order and limit order traders, reveals 
widening of the bid-offer spread when the flux of market order traders exceeds that of 
limit order traders. The variation of the spread with the fraction of market-order traders 
is investigated. The formula for asymmetric fluxes is applied to the two-liquid model and 
its predictions are found to agree with experiment. The critical point is approximately 
determined, and the ratio of the midmarkets for / = 0.0 and / = 0.5 (where / is the 
fraction of market-order traders) is calculated. 
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1 Introduction 



Financial markets are the object of considerable interest to an increasing number of physi- 
cists, because of their complex nature and the applicability of stochastic techniques to 
such dynamic systems. This is a relatively new area of interdisciplinary physics research, 
and is catching the attention of quite a few theoretical and condensed matter physicists, 
especially those working in the fields of complexity and chaos. The relationship between 
physics and finance is not something that started in recent years. In fact, as far back as 
1900, Louis Bachelier proposed the random walk model of the stock market. Historically, 
much of the initial work on Brownian motion and fractals originated from studies of stock 
market behaviour, before being absorbed by mainstream physics. In the 1960s and 1970s 
these ideas became very popular and eventually led to the famous Black-Scholes option 
pricing formula. 

Recently, many papers have appeared in which markets are treated as far-from- 
equilibrium dynamical systems. Econophysics is a quickly growing field of research, and 
work is being done on the scaling behaviour of exchange rates, "log-periodic" oscillations 
as crash precursors, dynamics of the interest rate curve, and market fluctuations, to name 
but a few examples. Some list of references (very incomplete) can be found in | 10[| . For 
more detailed discussion on physical approaches to financial markets and economy we 
refer the reader to recent publications fig] , [gj, Jll| which provide a lot of references. 

In this paper we present results of numerical simulations of two models - Minimal and 
Two-liquid, which had been formulated by two of us (D.E. and U.K.) in the paper on 
market microstructure ||10j| . In our approach the market microstructure is described by a 
diffusion-annihilation model. Different states of market are described by different possible 
state of a system in which diffusion controlled annihilation reaction takes place and the 
simplest one is the steady state. 

The first application of diffusion-controlled annihilation to financial markets was sug- 
gested by Bak, Paczuski and Shubik|| to describe market evolution. These authors in- 
troduced a series of models based on diffusion-controlled annihilation as a route towards 
recovering the observed Levy-Pareto "fat-tail" distributions which are said to describe the 
medium term evolution of the stock market. In BPS model market microstructure, i.e. the 
structure of reaction zone was not considred at all. For discussions about similarities and 
differences between BPS and our approaches see jl(J. Let us also note that in this paper 
we do not give references on economical and financial papers on market microstructure. 
The list of more than half a hundred important publications is given in [10|. 

In our model, we shall be studying market microstructure measured over short time 
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scales, so that sociological interactions between traders may be ignored. In the next few 



sections, we will describe the details of the models (§ |3.1| and § |3.2|) . 

Much research has been done in steady-state diffusion- driven annihilation reactions, 
especially in physics and chemistry. In applications in physics and chemistry, such re- 
actions are usually in three dimensions (e.g. in a solid or a liquid), but in our model of 
financial markets, we deal with one dimension only. In the initial paper by Galfi and 
Raczfl^], the properties of the reaction front in a system with segregated initial condi- 



tions were studied in the mean field approximation. The mean field approximation, which 
ignores higher, non-Gaussian corrections, known as fluctuations, works well for dynamics 
in three dimension, but breaks down in lower dimensions due to the important role of 
microscopic density fluctuations in one and two dimensions. Numerical simulations were 
performed in || |], |[ Analytical calculations by Cardy et. al. |LSJ [14], f| confirmed 

these numerical results. 

This paper is organised as follows. In the next section brief introduction into diffusion 
driven annihilation reactions (sometimes called DCR, i.e. diffusion controlled reactions) is 
given. Then applications of DCRs to market microstructure are discussed (mostly based 
on ||10|| ) in the section 3. In the last two sections the results of numerical simulations 
for Minimal and Two-Liquid models are presented. Some of this results confirm earlier 
known analytical or numerical results, but most of them are new - for example results 
about bid-offer spread, time to midmarket sale in Minimal model and critical ratio of 
market order traders in the Two-Liquid model. In the conclusion we discuss obtained 
results and new interetsing problems worth studying in a future. 



2 Diffusion driven annihilation reactions 

2.1 Basic principles of diffusion controlled reactions 

Diffusion controlled reactions, or DCRs, are reactions which include as one of the stages 
the transport of components which can be described by diffusion equations (see for exam- 
ple |11|). In DCRs, the rate of transport of reacting particles to within reacting distance 
of one another is the rate-determining step, and the rate of actual chemical combination 
is assumed to be much faster than that of transport by diffusion. This is a common 
situation that occurs not only in chemical reactions between atoms, molecules and ions, 
but also in cases involving electrons, lattice defects, quasi-particles, dislocations, etc. 

In classical chemical kinetics, the rate of chemical combination is assumed to be suffi- 
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ciently low compared with the rate of transport, that the uniform particle distribution is 
not disturbed and the system of two species may be considered to be in quasi-equilibrium. 
The rate of reaction is therefore a function of the state of the system, determined by 
the temperature, pressure and particle concentration. Diffusion space relaxation is fast 
enough so that diffusion is not the rate-determining step. In the classical approximation, 
therefore, we may characterize the dynamics by simple functions of state. 

This is not the case in the diffusion-controlled regime. Inhomogeneities arise because 
of local depletion of reactants in regions of high reaction rate. The overall rate of reaction 
is determined not only by the mean reactant concentration and other functions of state 
but also by the relative positions and velocities of the particles. The inhomogeneities 
have linear dimensions of the order of the inter-particle separation. In such cases, the 
notion of a local concentration is meaningless and the system must be described by a 
multi-particle distribution function. The search for the probability of particle collisions 
is much complicated by this. The mutual correlation in the positions of the particles 
is determined by the chemical reaction, making the probability of collision between two 
particles a complex function of the coordinator of all the particles in the system. In short, 
we have come up against the many body problem. Methods of quantum mechanics and 
statistical mechanics are widely applied in the search for a solution. 

DCR theory, which seeks a quantitative understanding of DCRs, consists of two stages. 
The first stage involves simple models based on single-particle distribution functions. The 
second stage incorporates the role of correlation effects. There are four basic assumptions: 

1. A System containing chemically active 'primary' reagents is in a state of thermal 
equilibrium. In a process consisting of a fast and a slow stage, the fast stage produces 
species which first achieve a thermal equilibrium, before these unstable species react 
further by the slow process to reach a chemical, and therefore thermodynamic, 
equilibrium. This is well illustrated by high energy radiation chemical processes, 
whereby radiation incident on a sample causes excitations to be produced in the 
molecules, forming ions, radicals and other reactive species. These species are non- 
uniformly distributed, usually in the form of tracks or cords. The system first comes 
to thermal equilibrium through exchange of heat energy between parts at different 
local temperatures. There then follows a progression towards chemical equilibrium, 
and thus thermodynamic equilibrium, through the ensuing chemical transformations 
of the unstable species. 

2. Initially, chemically active particles are assumed to be distributed inhomogeneously 
throughout the volume in the general case. 
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3. The main postulate of DCR theory is that active particles diffuse according to Fick's 
law, J = -DVn(r) where J is the diffusion current, D is the diffusion coefficient 
and n(r) is the local active particle concentration. This is applicable in the presence 
of strong macroscopic concentration gradients. However, in a spatially inhomoge- 
neous distribution, the notion of concentration is meaningless, as the characteristic 
length scales involved are of the order of the intermolecular separation. Here, we 
have to resort to analyzing diffusive mobility on a microscopic molecular level. We 
can imagine the particles situated on an imaginary lattice and hopping from site to 
site within a short time interval. After each jump, the particle comes into thermal 
equilibrium with the medium, from assumption f . The jumps are statistically in- 
dependent, represented by a random walk. This leads to the use of the theory of 
stochastic processes. 

4. Chemical reactions determine boundary conditions of the differential equations. The 
system is initially in a state of quasi-equilibrium. The chemical reactions do not 
change the inhomogeneous state of the system. The problem of how to take into 
account the limiting transport (diffusive) stages of the chemical transformation was 
solved by Smoluchowski by considering the probability for a certain particle to 
remain unreacted by a certain time. The dynamics of unreacted particles are still 
determined by the diffusion equations. Thus, the probability of reaction can be 
determined by the flux of mobile particles through the boundary assigned by the 
selected particle surface. Taking into account chemical transformations in DCR 
kinetics is reduced to the problem of applying boundary conditions to the diffusion 
equations. 

2.2 Applications of diffusion driven annihilation reactions 

Diffusion driven annihilation reactions have been well studied in physics and chemistry, 
partly because of their wide applications across many fields of these subjects. The ap- 
plications in chemistry are obvious; after all, chemistry is about the study of chemical 
reactions, a large number of which are diffusion-limited. There are also applications in 
physics, especially in the field of condensed matter where large numbers of particles are 
being considered. A newer, and perhaps less obvious, application is in economics, where 
the market microstructure determined by the short term behaviour of traders may be 
viewed as a form of diffusion through price space. We shall not go into it in any more de- 
tail here as the entire §[3] is devoted to this application of DCRs. Here are three examples 
of applications in chemistry and physics: 
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1. Atomic species A and B diffuse through a given volume with diffusion coefficients 
Da and Db- When A and B are a distance R from each other, they recombine in 
time dt with a priori probability -^rdt. Thus 

A + B -> C 

(inert) 

in a unidirectional reaction. If we have a bidirectional reaction where the product 
is not inert, 

A + B # C 

the reverse process has a decomposition rate of ^ye _/3C/ where /3 = 1/kT. Normally, 
U > but if the reaction is mainly that of decomposition, then U < 0. 

2. Fluorescence of certain liquid or amorphous semiconductors. The electrons and 
holes within the semiconductor diffuse separately and independently. When they 
are in close proximity to one another, they have a finite probability to recombine 
and emit gamma rays. Here, the unidirectional model applies. It also applies to the 
diffusion and decay of excitons at recombination ('scavenger') sites, an alternative 
and important mechanism for delayed illumination which involves only a single 
species. 

3. The structure of imperfect solids contains an excess of vacancies (missing atoms) 
such as those created by radiation damage. Vacancies diffuse throughout the solid 
until they recombine with interstitial atoms, or diffuse to the surface and effectively 
disappear. Alternatively, vacancies or interstitials created at the surface can diffuse 
to the interior. Diffusion-limited aggregation, whereby certain solids grow from 
vapour or liquid, has proved a popular subject of research. 

2.3 Uses of quantum field theory in diffusion-limited reactions 

Over the past two decades or so, it has been found possible, by many researchers, to use 
quantum field theory [|T7|j , || to study the motion, diffusion, recombination, and other 
dynamics of many body systems of non quantum mechanical objects, i.e. objects for 
which ApAx ~ H plays an insignificant role. It is used as a counting device, and is 
familiar to many physicists, especially particle and theoretical physicists. In quantum 
field theory, particles are constrained to execute random walks on the vertices of a space 
lattice in d dimensions. Quantum field theory allows a formal solution of the 'master 
equation' that governs many-body probabilities. To illustrate the method, we describe its 
application with an example. 
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The gist of the method may be demonstrated by considering the simplest example of 
all, radioactive decay. This is a zero dimensional problem, so we can proceed without 
worrying about diffusion operators. In a group of N radioactive atoms, we expect the 
number remaining undecayed at time t to be n(t) = N exp(— t/r) where r is the mean 
lifetime. n(t) is the expected number remaining undecayed, averaged over a large ensemble 
of such N atoms. It is known that initial conditions |PT|, p981] affect the short and medium 
term behaviour of the system, but the asymptotic behaviour for t 3> r is independent of 
initial conditions (as borne out by simulations in [§]). 

If 1/r is the rate at which a single nucleus decays, then, for the set of n undecayed 
nuclei, the rate at which a single decay occurs is n/r. The master equation for the 
probabilities is a linear differential-difference equation, 

dP{n\t) = {{n + l)P(n + l\t)- nP(n\t)}—, (2.1) 

r 

which incorporates the two unidirectional processes: decay into the state of occupancy 
n from n + 1, and decay out of it into n — 1. This master equation may be solved by 
indirection, by associating the state of n particles with the nth excited state of a harmonic 
oscillator pO| . Consider the harmonic-oscillator raising operator a* and its conjugate op- 
erator a. The commutation relation [a, a*] = 1 holds. If we define |0) as the ground state 
with zero occupation number, such that a\0) =0, the following relations may be obtained. 
We define \n) by 



This implies 



a*\n) = \n + 1), a\n) = n\n — 1). 

\n) = (oTIO), 
^a(a*) n \0) = ^(aT^IO), 
a*a(a*) n \0) = n(a*) n \0). 

It follows that a*a is the number operator. We now introduce the crucial concept of a 
probability state vector in which to embed the P(n\t): 

oo oo 

|*(f)) = £ P(n\t)\n) = £ P(n\t)(a*r\0). (2.2) 

n=0 n=0 

The initial condition is that there are iV atoms at t — 0, so P(n\0) = S n ,N, which is 
equivalent to the probability state |^(0)) = (a*) N \0). The right-hand basis states in this 
vector space (Fock space) are the (a*) n |0). Given 

l/nl(0\(a) n (a*) n '\0) = 8n,n> 
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they form a complete orthogonal set. However, the states are not normalized in the 
conventional manner. Instead, they are normalized with respect to a "reference" state 
(<S| = (0|e a which is merely a special case of a Glauber state (aS\ = (0|e aa (itself an 
eigenvector of the operator a* with eigenvalue a). 

The norm of any right-hand state |$) is defined in terms of (S\ by the inner product 
(5|$), so that a state is normalized if it satisfies («S|$) = 1. Given that (S\a* n = (S\ for 
all n > 0, and («S|0) = 1, each of the right-hand basis states (a*) n \0) is normalized in this 
fashion, i.e. (5|(a*) n |0) = 1 for any n > 0. Note that the probabilities remain normalized 
at all times t if they were normalized initially, provided the equations of motion conserve 
probability (which they always do): 

(SMt)) =Y. P (n\t)(S\(aT\0) = J2 P (n\t) = 1. 

n n 

Expectation values may be computed as follows. To evaluate the average number remain- 
ing at any time t > 0, we make use of the number operator a*a, 

(S\a*a\V(t)) = J2nP(n\t) = (n). 

n 

In general, the expectation value of a function of the number operator F(a*a) may be 
computed thus: 

(S\F(a*a)Mt)) = £F(n)P(n|f) = (F)(t). 

n 

The state vector |^/(t)) has all the properties of a generating function. First, it satisfies 
an elementary differential equation equivalent to Eq. (|2.1|) . Second, it yields individual 
P(n\t) by projection onto l/n\(0\a n . Third, it can be used to obtain the various moments, 
through contractions with (S\. The differential equation satisfied by is 

dt\v(t)) = — nmt)), (2.3) 

T 

in which Q = a*a — a, the dimensionless rate operator for this process, is sometimes 
referred to as the "quantum Hamiltonian" of the model. We illustrate how Eq. (|2.3| ) may 
be derived from the state vector equation fl2.2|) . The technique is to differentiate Eq. (|2.2j ) 
with respect to time, substitute in the master equation ( |2.1| ) and then try to arrange it 
into the form of a Schrodinger equation, as in ( |2.3|) : 



dP(n\t) 



\n) 



n=0 dt 

oo j 

YUn + l)P(n + lit) - nP(n\t)}-\n) (2.4) 

n=0 1 
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Now, we observe that 

oo 

(a*a-a)\ty) = (a*a - a) ^ P(n\t)\n) 



n=0 

oo oo 

nP{n\t)\n) - ^ nP{n\t)\n - 1) 

n=0 n=l 

oo oo 

nP{n\t)\n) - J2( n + l ) P ( n + MW 71 ) ( 2 - 5 ) 

n=0 n=0 



which we can identify with the right hand side of Eq. ( |2.4| ). Therefore, we can transform 
the master equation involving probabilities to a Schrodinger-like equation involving \^): 

dtim) = - {a * a ~ a) \v(t)) 

T 

T 

giving us Eq. (|2.3|). The solution is simply 

\V(t)) = e-^ rn \^(0)) = e~^ rn (a*) N \0), 

where the right hand side represents the initial condition of having precisely TV undecayed 
nuclei. The exponential in the solution may be expanded in series form and the whole 
expression substituted back into Eq. ( |2.3|) to check that it satisfies the differential equation. 

The above illustration of the method for dimensions can be extended to 1, 2 and 3 
dimensions, with the introduction of diffusion and annihilation operators. Although the 
example shown above is trivial, quantum field theory is a powerful method for solving for 
the dynamics of particles diffusing in a lattice, and allows analytic solutions to be obtained. 
It is beyond the scope of this report to discuss further examples of such applications, 



though the interested reader is referred to [17], [n] for a detailed review. 



3 The market microstructure of the interdealer bro- 
ker markets as a reaction zone of one-dimensional 
diffusion driven annihilation reaction 

3.1 The minimal model 

Because we are only interested in scaling properties, and not in modelling very detailed 
behaviour, the model is limited to only a few necessary parameters. This allows us to 
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Figure 3.1: A typical distribution of buyers (balls with letter B) and sellers (balls with 
letter S). Arrows describe "hops". 




B + S^O 



Figure 3.2: Annihilation B + S — > 



focus attention on the main model feature under examination. In ||10|| , it was proposed 



that the feature dominating double auction market microstructure was the great press of 
numbers of traders, so that the behaviour of, e.g. moments of measures of liquidity, were 
largely determined by the statistics of the traders. The required model is one in which a 
large number of individual traders interact according to any plausible dynamcics for the 
traders, and it is natural to choose the simplest one. Trader dynamics require that when 
a buyer and seller arrive at the same point in price space, they have agreed on a price, 
and must do a deal and vanish from the price space. And the simplest plausible model of 
trader motion is a random walk. Since we only wish to model markets in a steady state, 
we give these traders no net drift. 

This model of microstructure of the interdealer broker market can be mapped onto 
a statistical field theory, in which particles diffuse and annihilate. Thus, we associate to 
each buyer his bid price, and each seller his offer price, and imagine the two types of 



trade prices moving around in one-dimensional price space (Fig. |3.1|) . This price space is 
a discrete lattice, with a finite size. The traders, which may be thought of as particles, 
hop from one lattice site to the next with a certain probability (or they may choose to 
stay where they are). In this model, each trader trades in a standard size. When a buyer 
and a seller meet up at the same point in price space, they annihilate, leaving the market 



(Figs. |3.2| and |3.3| ). The reaction may therefore be represented by B + S — > 0, which is a 
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Figure 3.3: A new distribution of buyers (balls with letter B) and sellers (balls with letter 
S) after an annihilation (trade) has occurred. Note that the bid-offer spread fluctuates. 
Here it is bigger than in Fig. |3.1] 



diffusion controlled (or diffusion driven) annihilation reaction. This type of reaction has 
been studied extensively in the physics and chemistry literature |jTB| . 

We treat the case of a quiescent market, one which is an approximately steady state 
market, and we define our evolving state measure, at any time t, to be that probability 
measure which is the steady state solution of our diffusion-annihilation dynamics. As 
a first approximation, we allow traders to enter only at the ends, as though the buyers 
start by bidding very low and the sellers offering very high. One can solve this model by 
numerical simulation. However, because of its similarity to many well-known analytically 
soluble models (for example, in fll7l), we can go much further. Approximate analytical 
methods are available from Cardy et alO, O. ffl. 



3.2 The two-liquid model 

The minimal model is a suitable approximation of the quiescent periods of a market's 
behaviour, but when it is in the process of a long crash, it does not look quiescent and 
the model breaks down. It is not sufficient merely to cause an imbalance in the rates of 
insertion of buyers and sellers, because such a method would not reproduce the well-known 
phenomenon of the widening of the bid-offer spread. It is believed that the missing element 
is market order traders. During a crash, most traders, desperate to sell, attempt to hit 
the best bid directly, as a market order. In such a scenario, although limit order trading 
continues as before, the large proportion of market order trades cannot be neglected and 
indeed they have a significant effect on market dynamics. In the two-liquid model, we 
introduce a second type of trader, the market order trader, who tries to hit the best 
bid/offer directly, but without displaying his/her bid/offer price on the trading screen. 
These market order traders trade with limit order traders, who display their prices on the 
trading screen, but not with one another. If limit order traders are represented by B and 
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S and market order traders by B' and S', then the reactions B + S — > 0, B 1 + S — > 0, 
and 5 + 5" — > can occur, but noi 5' + 5" — > 0, since market order traders are invisible 
to one another. This model is able to reproduce the familiar widening of the bid-offer 
spread during a crash. 



3.3 Dimensional analysis 

Because we have relatively few degrees of freedom (parameters) in our model, we can use 
dimensional analysis to obtain the scaling forms of many market quantities, as a check on 
intuition. Starting from the operator evolution equation, we may differentiate it to yield an 
operator differential equation. This partial differential equation consists of one coefficient 
only, namely the diffusion coefficient D, with dimensions x 2 /t, i.e. D ~ (dollar) 2 /sec (see 
§|A.1| for derivation of the diffusion coefficient). In addition, there is another dimensionful 



quantity in the boundary conditions, J, the rate of insertion of traders at the boundaries, 
which is equal to the deal rate in a steady state market, with dimensions J ~ l/(sec). 
From D and J we can construct quantities with any dimensions containing (dollar) and 



(sec), since (dollar) ~ yD/J and (sec) ~ 1/J. Therefore, the expectation value (X) 
of a general quantity X with dimensions [X] = (dollar) m /sec n must be proportional to 

(D I ' J) m / 2 J n ~ ]jrn,/2 jn-m/2 

So far, we have ignored, in our dimensional analysis, of the lengths scales L, the size of 
price space, and 5S, the lattice spacing. These two lengths represent the upper and lower 
length scales in our model. We are justified in doing so by the fact that if we imagine a 
Fourier analysis of the dynamics, most of the 'action' would be taking place away from 
these extreme length scales, in the intermediate region between the two. Ultimately, 
the justification comes from numerical simulations. However, these length scales do add 



logarithmic corrections to the scaling laws, as shown by Cardy et al|15], [LJ, |4], [L(|. This 



is described in 53.6 



3.4 Market parameters 

The model allows us to calculate the relationships described in the introduction, but also 
many more. In general, whatever initial conditions we start with, in the asymptotic limit, 
we will reach a steady-state, with a reaction front, where the buyers meet the sellers. At 
the centre of the reaction front is the best bid, best offer and midmarket. Beyond the 
best bid/best offer, we expect to see an increasing density of buyers/sellers. Traders are 
inserted at the edges at a rate equal to the deal rate. 
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The statistics describing the state of the system can be read off a trading screen. 
The best bid locates the top of the lower edge and the best offer the bottom of the 



upper edge (see Figs. |3J] and p.3| ). Although there is a large number of interesting 
parameters, in this project we shall be restricting our attention to a small subset of them: 
best bid B(t), best offer 0(t), bid-offer spread Spr(t) = 0(t) — B(t), and midmarket 
M(t) = (1/2) (B (t) + 0{t)). Note that B(t), 0(t), Spr(t) and M(t) are not Markov 
random variables, and do not satisfy a stochastic differential equation, because they are 
subject to jumps. In addition to the above parameters, we shall also be interested in 
recording the instantaneous deal rate DR and the total number of traders in the market 
NUM. 



3.5 The scaling laws 

Financial markets, like many other systems, exhibit scaling laws which are clearly observed 
by practitioners in the field. These scaling laws are related to the concept of liquidity, 
which is measured in several different ways. It may be measured by the deal rate J, 
the time to midmarket sale, the bid-offer spread, and the trader densities near the best 
bid/offer. Here we shall concentrate on several scaling laws which were introduced in [|K] 



Bid Offer Spread 

Spr = (0(t) - B{t)) (3.1) 
Trade price variance (or midmarket variance) 

w 2 = (M 2 (t)} - (M(t)) 2 (3.2) 



• Fluctuations in the Bid-Offer Spread 

(ASpr) 2 = vav(Spr) = ((0(t) - B(t)) 2 } - (0(t) - B(t)} 2 (3.3) 

• Instantaneous Dealing Time, and its Fluctuations 

We can define Ts as the time between successive trades. If after a time, more than 
one trade occurs in the same time step, then, in order to have a valid continuum limit 
for the quantity, we must consider the time between the multiple trades themselves 
to be zero. We may visualize it thus, assuming time is continuous: one trade occurs 
at t — 0, and then no trades occur for some time. The next trade occurs at t — 9, 
and the next at t = 9 + 5. These are distinct times, if we use continuous time, and 
we would say that the two values of t$ are 9 and 5. If we now allow 5 to tend to 
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zero, we have 9 and as the values of Ts, but we see that such a configuration is 
equivalent to the discrete model of time, where we have no trade for 9 time steps 
followed by 2 at the same time. Thus, the correct way to deal with "simultaneous" 
annihilations is to treat them as separate ones, separated by zero time. 

We can define its fluctuation as 

Ar s 2 = var(r 5 ) = (r s 2 ) - (r s ) 2 (3.4) 

Asymmetric fluxes 

It is possible to consider the minimal model with different fluxes for the buyers 



and the sellers. In [|T_0j the following theoretical result, valid for "small" times, was 
found: 

-2DAJ . . 

x (t) jj—t, (3.5) 

where xo(t) is the position of the midmarket, which varies with time t. The flux of 
sellers from the upper boundary is J + A J whilst the corresponding flux of buyers 
from the lower boundary is J — A J. 

We expect Spr(J) and w 2 to tend to as J — > oo. From dimensional analysis, Spr and 
w have dimensions of dollars, so Spr(J) ~ J[D/J), w(J) ~ J(D/J) consistent with 
intuition. Using the results of [JTR [TJ|, ffl, we see that there is a logarithmic correction 
for w. It is also possible to measure fluctuations in the spread, ASpr or var (Spr). This 
quantity should have the same scaling law as Spr, since it has the same dimensions. Time 
to midmarket sale t$ is expected to go as ~ 1/J, from dimensional analysis; in fact, its 
average should be equal to 1/J. Similarly, its fluctuation, Ats should have the same 
scaling law. There are many other scaling laws which we will not go into here as they will 
not be tested later; for more information, see |10| . 

3.6 Analytical results 

Although numerical results are more precise, it is desirable to work out some scaling laws 
for the model using analytical means, as these methods supplement numerical simulations 
by providing intuition. Cardy et. al.JTR 113, £1 has developed an approximation scheme 



known as mean field theory, replacing the evolution equations for the operators with a 
partial differential equation for the density configuration with the greatest probability 
mass, ignoring the "fluctuations" . There is no system of higher order corrections or error 
estimate with this method. In spite of its limitations, its predictions have been found to 
coincide with numerical simulations |L5], [14], |], [L| ^, ||. It is beyond the scope of this 
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paper to describe in detail the derivation of the analytical results, so we shall simply state 
the most important result for our simulation — the logarithmic correction to the scaling 
law for to: 

_ r in(cL/w) l 1/2 
W ~ [ <J/D) _ 

This was derived by the addition of an effective noise term to the differential equation 
satisfied by the density difference of buyers and sellers. However the numerical value 
of constant c is unknown from analytical calculations and we shall find it later from 
numerical data. 

The second analytic result we shall be using relates to the speed of the moving mid- 
market in the case of asymmetric fluxes of buyers and sellers. The result here is (for more 



details see page 34 in |[T0| ) 



x (t) 



-2DAJ 



JL 



-t 



and was derived using a time-dependent non-stochastic part in the density difference 
equation. It is only valid for small times, such that T < L 2 /D. 



4 Computer simulations of the minimal model 

4.1 Description of simulation 

We shall start from a Monte-Carlo simulation of the minimal model. 

The simulation starts by randomly inserting traders into either half of price space to 
set up a random initial configuration. The main loop is entered. Traders are inserted at 
the edges, and then each type of trader (buyer and seller) is allowed to diffuse separately 
and annihilate. Each trader hops to the left or the right at every time step with probability 
0.5 (p = 1). This leads to a diffusion coefficient of 1/2, because D = a 2 /2r (see §A.1|) , 
and a = r = 1 in our units. The annihilation routine is called after each diffusion 
routine to ensure that pairs of mutually annihilating traders do not cross over each other. 



This is in accordance with the procedure outlined in the appendix of JT0[, where the 
time-evolution operator is constructed from diffusion and annihilation operators, with 
annihilation occurring after each diffusion of buyers and sellers (cf. [|| which has a slightly 
different way of doing it). The tasks performed within each time quantum r may be 
summarized as: 

1. Diffuse buyers 
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2. Annihilate overlapping traders 

3. Diffuse sellers 

4. Annihilate overlapping traders 

5. Call getstatsO routine to store market parameters 



This method we shall call MINIMAL. At the end of each time step, the getstatsO 
routine is called, which obtains interesting statistics from the state variable, such as 
best bid, best offer, bid-offer spread, midmarket, instantaneous deal rate, etc. Thus the 
diffusion-annihilation process is repeated over and over again. Now, not all the data 
generated at each time step is recorded, simply because it would be impractical and 
unwieldy to have a set of two million results to deal with. Instead, results such as bid- 
offer spreads and midmarkets are averaged over a large number DISP_INT of time steps. 
In the end, we would like to take the average^ and variance of all the results, not just the 
results averaged over DISP_INT steps. At first, one might be inclined to think that the 
process of averaging over some steps causes one to lose information. In fact, it is possible, 
by simultaneously keeping track of the mean of the squares of the quantities, to allow the 
mean and variance of the entire ensemble to be recovered, as though we had kept track 
of every single result (for proof see Appendix § A..2 ). 



The data produced from the program were recorded in a data file. Note that in the 
simulation, the quantities considered were dimensionless, and may be defined as follows: 

T = Aqr 
L = N 2 a 
Jr = N 3 

where r is the time between each step in the simulation (the time quantum), T is the 
total time of the simulation, L is the width of price space, a is the lattice spacing of 
price space, and J is the rate of insertion of traders (no. of traders per second). N\, N 2 , 
and N 3 may be thought of as the dimensionless counterparts of T, L and J, respectively. 



§Note that the average taken of the results generated by computer simulation is a time average, whilst 
the "average" used in the scaling laws and the analytical results ((. . .)) is an ensemble average. These two 
are generally not the same. However, for a stationary ensemble (i.e. one that has no preferred origin in 
time), the ergodic assumption, that the system will in the course of a sufficiently long time pass through 
all the states accessible to it, leads to the fact that the time average is equal to the ensemble average. 
The financial market, when in equilibrium, is well-approximated by a stationary ensemble, and so the 
equality of the ensemble and time averages applies. That is why it is acceptable to compare theoretical 
ensemble averages with the time averages obtained from numerical simulations, which are the same as 
ensemble averages. See [M for proof. 
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We may simply choose our length and time units such that a = r = 1, so that lengths 
are measured in units of the length quantum a and time intervals in units of the time 
quantum r. In these units, Ni = T, N 2 = L and N 3 = J, and so we may revert to the 
natural units, which have been made dimensionless by a judicious choice of units. 

4.2 Preliminary results for L = 200 using MINIMAL 

The simulation of the minimal model was initially performed for 100,000 time steps, with 
averaged results recorded every 100 or 1000 steps. The number of steps in the simulation 
was limited by the length of time taken to complete it: for each run of 100,000 steps and 
one particular value of J, the run-time would be between 1 and 2 hours. The time taken 
to obtain an entire set of results was a few days. There were 200 points in price space 
(L = 200), and the initial condition consisted of 100 buyers and 100 sellers distributed 
randomly on either side of price space, with a gap of 20 empty price points in the middle. 
^From this condition, the system was allowed to evolve according to the diffusion and 
annihilation laws outlined earlier. Every DISP_INT time steps, the following averaged 
quantities were recorded: best bid (B), best offer (O), best bid size (BS), best offer 
size (OS), midmarket (M), bid-offer spread (Spr), deal rate (DR), and total number of 
traders in the market-place (NUM). 

It was found that the system came to equilibrium after about 50,000 time steps. The 
total number of traders in the market (NUM) was used as a suitable measure of the state 
of equilibrium of the system, because it was discovered that NUM grew with time and 
tended asymptotically to an equilibrium value (Fig. [41]) . We have also included a picture 
of how the bid, offer and midmarket vary with time in the simulation (Fig. |4.2j ). 

Thus, averages of quantities were only taken over the range of time steps for which 
the system was in equilibrium, in this case, the last 50,000. 

The scaling laws we are interested in are the following: 



To this end, we plot graphs of In w 2 against In J, and In Spr against In J. We can see that 
each run of the simulation produces one point on a graph, since we average over a long 
time to obtain mean values for our market parameters. To obtain any reasonable graph 
requires at least 10 points or so, hence the long times taken to obtain results. 




(4.1) 



(4.2) 
(4.3) 
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Figure 4.1: Graph of NUM against time for J = 4 
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Figure 4.2: Variation of bid, offer and midmarket with time for J = 1 
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Figure 4.3: Graph of \nw 2 against In J 
for J > 1 and L = 200 




Figure 4.4: Graph of In Spr against In J 
for J > 1 and L = 200 




Figure 4.5: Graph of lnvar(S'pr) 
against In J for J > 1 and L = 200 




Figure 4.6: Graph of Jw 2 against hiw 
for J = 1 and L = 200 



The results obtained (Figures [4.3| , [4.4| , [4.5| and |4.6| ) were of a poor quality and did not 
yield the expected straight lines. The values of J used were 1, 2, 4, 8, 16, 32 and 64. One 
can see that these values of J are too high and that in this regime, the market is very 
dense and effectively frozen. This is not the regime described by the minimal model, and 
it is not surprising that the results do not agree with the predictions. From the data, 
one can see that w < a, var(Spr) < a, and Spr ~ a. The reason for the poor results is 
the following: the scaling laws, derived using dimensional analysis, relied on the fact that 
the length scales of meaningful quantities were far from L and a. Thus, for dimensional 
analysis to hold, and the scaling laws obtained thereby to be valid, we require that 

d«w, Spr, var(Spr) L. (4.4) 

This means that J must be much less than 1, since all the length quantities are of order 
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y[{D/J). 



4.3 Successive improvements to preliminary results 



We outline a sequence of improvements to the simulations as follows: 



Method L T 



Range of J 



MINIMAL 200 1 x 10 5 

MINIMAL 200 1 x 10 5 

MINIMAL 200 1 x 10 6 

MINIMAL 1000 2 x 10 6 



J > 1 

J < 1 

J < 1 

J < 1 



The average time taken to obtain each set of results was about 3 or 4 days. The simulations 
were repeated for fractional J, from 0.001 to 0.1, for L = 200. The results obtained 
(MINIMAL, L = 200, T = 1 x 10 5 ) were better than before, and the points lay more 
closely on a straight line. As the simulation time is proportional to J (see § |4.13|) , it 
becomes possible, with J reduced by several orders of magnitude, to increase the number 
of time steps per run to 1 million whilst keeping the total duration reasonable, with 
the statistics being taken over the last 950,000 steps. Thus, it was decided to do the 
simulations again, but with 1 million steps. The results obtained (MINIMAL, L = 200, 
T = 1 x 10 6 ) were even better, though the graphs were still not very convincing, especially 
the one of Jw 2 against \nw, which should have been a straight line, according to Eq. (|4.1| ), 
but in fact the points were so scattered that it was difficult to draw any definite conclusion. 

The smallness of L (200) restricted the choice of J too much: it was desired to keep w 
at least a factor of 5 away from either a or L, so the range of allowed w was from 5 to 40, 
which was just one order of magnitude. Thus, it was resolved to increase L to 1000, which 
expanded the range of permitted w to 5 < w < 200. However, the increase in the size of 
the system caused a corresponding increase in the time taken for the system (especially 
trader number NUM) to reach equilibrium. It was noticed that quasi-equilibrium was not 
attained until after about 1 million steps. The duration of the simulation was therefore 
increased to 2 million steps, with only the last 1 million steps used for obtaining statistics. 

4.4 Encouraging results using MINIMAL for L = 1000 and T = 



Many different values of J were used, all of which produced values of w and Spr roughly 
within the range specified above, with most interesting lengths being at least a factor of 5 
away from either a or L. Here, the initial condition consisted of 50 buyers and 50 sellers, 



2 x 10 
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Figure 4.7: Graph of lnw 2 against In J for L = 1000 



separated by a gap of 20. One of the things learned from these simulations (obtainable 
by intuition, and supported by ||) is that the asymptotic behaviour of the market does 
not depend on the initial conditions. 

38 different values of J were used, producing a set of 38 points on the graphs (Figs. [0|, 
4~8| , [4.9| and 4.10j) . It was noticed that for the first three graphs at least, the points fit 
a line better towards the bottom right of the graph, while the points towards the top left 
tended to be more scattered. It was likely that the points towards the top left (small J 
and large w 2 ) were less accurate because of the possibility of the reaction front crashing 
into either edge of price space. In addition, for graphs of lnw 2 , In Spr and lnvar(Spr) 



against In J, the analysis in § [4.5| requires that lnw IncL be satisfied, and this is not 
satisfied for the points towards the top left, for which In w ~ 5. Therefore, we must choose 
a small section of the bottom right of the graph and obtain the gradient and intercept 
from that. 

Thus, it was decided to take only the last 20 points and fit a line through them, in 
the hope of acquiring more accurate statistics (Figs. |4.11| , [4.12| , [4.13| and |4.14j) . Here is 
a summary of the results, for the first set of data (all 38 points): 
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Figure 4.8: Graph of In Spr 2 against In J for L = 1000 




hJ 



Figure 4.9: Graph of Invar (Spr) against In J for L = 1000 
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Figure 4.10: Graph of Jw 2 against haw for L = 1000 
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Figure 4.11: Graph of Into 2 against In J for L = 1000 (last 20 points) 
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Figure 4.12: Graph of In Spr against In J for L = 1000 (last 20 points) 



b- 4 

to * 

i 

"13 

I 3 



2 
1 




R ■ square = 1 ttpts = 20 
y = -Q744 + -1x 





N 



-8 -7 



-5 -4 

. . In.J. . . 



-1 



Figure 4.13: Graph of lnvar(Spr) against In J for L = 1000 (last 20 points) 
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Figure 4.14: Graph of Jw 2 against \nw for L = 1000 (last 20 points) 



Figure Abscissa Ordinate 



Gradient 



Intercept 



R 2 



4.7 



4.8 



mi 



4.10 



In J 
In J 
In J 

lnu> 



lnw 2 
In Spr 2 
In var (Spr) 
Jw 2 



-0.814 ±0.015 
-0.925 ±0.014 
-0.924 ±0.014 
-0.174 ±0.013 



0.41 ±0.26 0.988 

0.903 ± 0.24 0.992 

-0.38 ±0.25 0.991 

1.006 ±0.089 0.844 



and here is the same for the second set of data (last 20 points): 

Figure Abscissa Ordinate Gradient Intercept R 2 



1.11 


In J 


In w 2 


-0.919 ±0.016 


0.052 ±0.1 


0.995 


4.12 


In J 


In Spr 2 


-0.9904 ±0.0024 


0.592 ±0.017 


0.99989 


4.13 


In J 


In var(Spr) 
Jw 2 


-1.001 ±0.003 


-0.744 ±0.023 


0.99980 


4.14 


In w 


-0.107 ±0.024 


0.888 ±0.076 


0.523 



It must be noted that for all the graphs, with the exception of Jw 2 against Inw, the 
correlation coefficients are all 0.98 or above, which means that the results strongly support 
a scaling relationship of some sort, the exact details of which we shall now discuss. 
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4.5 The w 2 scaling law 



Now that we have some reasonably detailed results, we may begin our analysis of them 
in earnest. The fact that the plot of Inw 2 against In J is a straight line with a gradient 
close to -1 suggests that w 2 is inversely proportional to J, which is correct according to 
Eq. (|4.1|) if we ignore the logarithmic correction. This is simply a consequence of the 



dimensional analysis in §[3]3 



Next, we seek to observe the logarithmic correction. Consider the plot of Inw 2 against 
In J. From Eq. ( [4.1D , and from D = a 2 /2r = 1/2, we can write 



In w z 



- In J + In 



2n \ w 



— In J + In ( — [ln(cL) — In w 
1 2n 



In J + In < — ln(cL) 
2ir 



\nw 
In cL 



- In J + In 



+ In 



In cL 



Here, we make the approximation ln(l — x) ~ — x, ignoring the higher terms —x 2 /2 — 
x 3 /3 — . . ., which is valid for x <C 1, on the third term on the right hand side. In this 
case, the approximation is valid if Inw <C In cL. Letting \nw = 1/2 Inw 2 , and grouping 
terms of Inw 2 on the left hand side, we obtain, 



1 + 



21ncL 



In w z 



- In J + In 



In cL 
2tt 



For 2\ncL ^> 1, we can make the further approximation of yl 
yielding 



i 



Inw 



2\ncL, 



-a In J + a In 



In J 



In cL 

2tt 



1 



1 



2 In cL 



In 



l 



2 In cL 



In cV 
2tt 



2 In cL ' 

(4.5) 
(4.6) 



which defines a. If we assume that c ~ 0(1) and take L = 1000, then the logarithm in 
a is positive, making a slightly smaller than unity. From measuring a, we may therefore 
deduce the value of c. In this way, it is possible to see the subtle effect of the logarithmic 
correction. Note that the second approximation is not strictly necessary and one can use 
the full expression for a if one wishes, and this is indeed what we shall do. 
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We make use of the graph of the last 20 points because we want the approximation 
In w <C In cL to hold, as otherwise the foregoing analysis would be invalid. From Fig. |4.11 , 
and using the full expression for a (i.e. without making the binomial expansion for the 
reciprocal) we find that the gradient of the graph is —0.919 ± 0.016, making a = 0.919 
and c = 0.29, but because c is inside a logarithm, the uncertainty in a is magnified in c, 
so that it lies in range 0.11 < c < 1.33. 

We can also obtain a value for c from the intercept, which is — 0.052±0.1, another wildly 
imprecise quantity. This method of calculating c gives 0.38, with limits 0.20 < c < 0.75. 
The intercept appears, therefore, to provide a more precise determination of c, though, 
as we shall see later, that is not always the case. 

There is another technique one can try in extracting information from the graphs. 
Instead of allowing only c to be determined from the graphs, we may attempt to determine 
the constant 1/tt from the gradient and the intercept, since each graph has two degrees 
of freedom. Therefore, we replace 1/tt with a variable, say, Q, and the equation becomes 

2 A 
\nw = —aln J + aln (—In cL I , 

and we may determine Q using the data. Here, we have Q ~ 0.333, with the uncertainty 
0.236 < Q < 0.453. Note that, the value predicted from theory is 1/tt = 0.3183, which 
lies well within the experimental error. This is quite encouraging as it agrees with the 
result derived theoretically (in spite of its 33% uncertainty!) 

It is possible, for the sake of comparison with § [4.6| and §4.7| , to analyze the data for 
w 2 as though its scaling law had no logarithmic correction. It is not really appropriate, 
because of the relatively pronounced (-0.919) deviation from -1, but we shall quickly do 
it for subsequent comparison. Assuming that w 2 is proportional to D / J with constant of 
coefficient A, 

2 AD 
W = ~ 

we can take the log of both sides 

In w 2 = In \D — In J 

and calculate A from the intercept, knowing that D = 1/2. The intercept is —0.052 ±0.1, 
giving us A = 1.90, within the limits 1.72 < A < 2.10, or 1.9 ± 0.2. 

We have another graph from which we may obtain information: the graph of Jw 2 
against \nw. For this graph, it is prudent to take the one for all 38 points, since no 
approximations have been made in the scaling law we are trying to verify, and the more 
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points we have on the graph, the more confident we may be of its statistics. It is simply 

Jw 2 = ^(lncL-lnw) (4.7) 



which comes from a trivial re-arrangement of Eq. (|4.1|) . We have once again allowed I /it 
be determined from the graph. The gradient is —0.174 ± 0.013, from which we may infer 
Q ~ 0.348 ± 0.026. This is close to the predicted 1/ir = 0.318, although it actually lies 
outside the range of experimental error (the discrepancy is 9%). We can, with reasonable 
confidence, now set Q to its theoretical value, 1/tt. From the intercept, 1.006 ± 0.089, 
we may infer c = 0.556 with the limits being 0.32 < c < 0.97. These limits are slightly 
different from those found from the first graph (Fig. |4.11|) , but not by much. Both methods 
give the same order of magnitude estimate c ~ 0.5. 



4.6 The Spr scaling law 



We shall make use of Figs. [4.8| and |4.12| , the latter consisting of the last 20 points of the 
former. It is interesting to notice that the correlation coefficients for both graphs are very 
high, and indeed, that of the second graph is so close to 1 that it was rounded up to 1 by 
the spreadsheet that produced the graph. Once again, the straight lines with gradients 
so close to -1 confirm the approximate Spr ~ J(D/J) law. It is not known (see |T0|j) 



whether there is a logarithmic correction to Spr or not. If there is, one might expect it 
to be of the form 



2 aD , / sL 



Spr = — — In 



J \Spr J 

in analogy with Eq. ( |4.1| ). We have left the constant in front of D/ J as a parameter (a) 
to be determined. 

To look for the logarithmic correction, we can perform a similar analysis as we did in 



§ |1.5| , making the approximation In Spr <C IncL, to obtain 

In Spr 2 = —a In J + a In ^— In sL^j 

where 

I = i + 1 



a 2 In sL 

To ensure the validity of the approximation, we shall use the graph of the last 20 points 
only (Fig. |4.12|) . There is another reason for putting more faith in the statistics yielded 



by the second graph: if one looks carefully along the points of the first graph, one sees 
that the points towards the lower right of the graph lie along a line that does not quite 
coincide with the line of best fit constructed from all of the points. This is because 
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of some deviation from linearity towards the top left of the graph. It would be wise, 
therefore, to take only the results towards the bottom right, as in Fig. [4.12| . The gradient 



is -0.9904 ± 0.0024, which would imply s - 2.5 x 10 19 , 7.6 x 10 14 < s < 8.8 x 10 26 . 
The uncertainty covers 12 orders of magnitude! The intercept is 0.592 ± 0.017, giving us 
a ~ 0.07, with limits 0.054 < a < 0.087. We shall try to understand the meaning of the 
values obtained. Consider a rearrangement of the scaling law 



Spr 2 = — — 

■J 



In s + In 



Spr 



The largest possible value for \n(L/Spr) occurs when Spr is at its smallest, around 5. 
Thus, the maximum value for the second logarithm is 5.3 (minimum is about 1.6). The 
first logarithm, however, is Ins, which is 44 (34 < Ins < 62). Therefore, the right hand 
side is dominated completely by the large constant In s, and the logarithm involving L 
has almost no effect on the scaling law. This is a clear sign that there is no logarithmic 
correction, since the variation of the second logarithm has negligible effect on the vari- 
ation of the sum of the two logarithms, on which the scaling law depends. This can be 
approximated by ignoring the variation of the second logarithm with Spr, and taking 3 
as the approximate average value for it, leaving In s ~ 47, with limits 37 < In s < 65. 

Spr 2 = Ins, 



which is, 

Spr 

where the coefficient \x = a Ins ~ 3.3. 

Now, if we assume that there is no logarithmic correction, we can try to extract from 
the data the constant of proportionality for the scaling law 

c 2 l*D 
bpr = ——. 

u 

Taking logs, we have 

In Spr 2 = ln(/iZ)) — In J. 

The gradient agrees well with the predicted value of — 1. We know that the diffusion 
coefficient D is 1/2, so from the intercept, \i ~ 3.615, with limits 3.55 < \i < 3.68, or 
3.615 ± 0.065. Thus, the two methods, one assuming a logarithmic correction, and the 
other not, produce results which agree with one another (the 3.3 obtained earlier has an 
associated uncertainty which is rather difficult to work out, so we will not do that here). 
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4.7 The vai(Spr) scaling law 



We look at Figs. 43 and |4.1ij . Both have extremely high correlations (almost perfect 



correlation for Fig. 4.13) , so we can have confidence in the results. The results are 



unambiguously linear, confirming the scaling relationship var (Spr) ~ y(-D/J). The only 
thing remaining is the determination of the existence of logarithmic correction. 

The gradient of the graph is —1.001 ± 0.003, which puts the exact value of -1 within 
the range of experimental error. Even if there is some deviation from -1, it is minute. 
This is strong evidence against the existence of logarithmic correction for vsa(Spr). There 
is no point in trying to apply the same analysis as before to find out the constant inside 
the logarithm, because here, the gradient can lie on either side of -1, causing the value 
of that constant to be anywhere between the exponential of a large negative number and 
the exponential of a large positive number, i.e. many orders of magnitude. The fact that 
the constant can assume these extreme values implies that the logarithm of the constant 
would dominate the scaling law, when it is either very small, in which case the logarithm 
of the constant would be a large negative number, or very large, in which case the log 
would be a large positive number. It is fair to conclude that Spr appears from the data 
not to have any logarithmic correction. 

We may proceed to find the constant of proportionality, z/, defined by 

var (Spr) = — — . 

'J 

Taking logs, 

In var(Spr) = ln(^D) — In J. 

The gradient has already been verified to be -1. The intercept is —0.744 ± 0.023. There- 
fore, v ~ 0.95, subject to 0.93 < v < 0.97, or 0.95 ± 0.02. 

4.8 Calculation of c for different values of L 



^From Eq. ( [4.6|) , we found that plotting graphs of In w 2 against In J allowed the logarith- 
mic parameter c to be determined. After further discussion, it was decided to investigate 
this parameter and try to see whether and how it varies with L. To this end, simulations 
were run for different values of L, within the correct regime for J, such that a <^ J ^ L. 
The observance of this regime was now especially important, owing to the ln(l — x) ~ —x 
approximation we made in our derivation of Eq. ( |1.6| ) . 

The results obtained were rather poor. Correlations were low and the value of c did 
not appear to show any kind of relationship to the value of L. We shall not present these 
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unhelpful results here. Instead, we shall describe a new method of simulation, which will 
be used to obtain better results. 

4.9 MINIM AL1: a new method 

If we take a detailed look at the variation of something like Spr with time, it will be seen 
that the MINIMAL method of simulation, which we have been using from the outset, 
does not allow Spr to go to zero. This is because of the sequence of steps undertaken at 
each time quantum r. Only after both species of traders have diffused and annihilated do 
we take measurements of the system using getstats(). In practice, annihilations in the 
marketplace are not instantaneous, and for some of the time at least, overlapping buyers 
and sellers can exist. This suggests an alternative method of simulation, MINIMAL1, 
based on the following 8 steps: 

1. Diffuse buyers 

2. Call getstatsO routine — overlapping traders possible 

3. Annihilate overlapping traders 

4. Call getstatsO routine — no overlapping traders 

5. Diffuse sellers 

6. Call getstatsO routine — overlapping traders possible 

7. Annihilate overlapping traders 

8. Call getstatsO routine — no overlapping traders 

This method represents a more detailed observation of the market as it includes interme- 
diate overlapped states which were ignored in MINIMAL. MINIMAL1 allows the bid-offer 
spread to vanish just before an annihilation, and can be said to produce a more accurate 
profile of the evolution of the market. However, equal weight is given to the statistics 
obtained at each of the even steps shown above. Whether this is a good picture of real-life 
trading is another matter altogether, which we shall not discuss here. 

4.10 MINIMAL1 in action 

The test of this model was to see whether it would produce better results than MINIMAL. 
We still used the method of averaging over DISP_INT time steps to obtain our statistics, as 
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Figure 4.15: Graph of bid, offer and midmarket with time, for L = 100 



two million time step results were too many to deal with easily. However, it was possible, 
by a further modification to the program, to produce code that instead of averaging over 
DISP_INT, wrote the state of the system at points 2, 4, 6 and 8 in the above sequence of 
steps, so that each getstatsO call was followed by the recording of the data in a file. 
This allowed one to produce a graph of how various quantities such as B, O and Spr 
varied with each and every time step. Besides being a satisfying thing to look at, such a 
graph gave a visual demonstration of the simulation at work, and an assurance that the 
simulation was working as it ought (Figs. [4. 15| and |4.16| ). The graphs show that bid and 
offer can meet at the same point at certain times, just before the annihilation of traders at 
the reaction front. Fig. [4. 1 6| demonstrates this especially well. Spr gradually gets closer 
to zero by the diffusion of the traders at the reaction front, until eventually it reaches zero 
and annihilation occurs, after which the best bid and best offer immediately 'snap back' 
to the position of the next best bid/offer. This is precisely the way the model is expected 
to work. Satisfied that our new simulation is functional, we press on and analyze results 
obtained thereby. 



4.11 MINIMAL 1 results 



We shall perform the same analysis for MINIMAL 1 results as we did for MINIMAL. To 
ensure the approximate validity of the approximation \nw <C IncL, we shall take only 
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Figure 4.16: Graph of Spr with time, for L = 100 



100000 



results for which Inw < 3, because if we take c ~ 0.5, then for L = 1000, IncL ~ 6. 
Values of Inw ~ 3 cannot really be considered small compared with IncL, but if we 
restrict our range any further, then we would be making use of a small number of results 
(10 points or so) in which case reliability of results would suffer. Thus, we will use points 
for which In it; < 3, bearing in mind that the points towards the top left of the graph are 
less accurate. 

It was decided that the last 12 points would be used to plot a graph from which 
the gradient and intercept would be obtained. Two factors influenced this decision — the 
desire to have as many points as possible to increase reliability of the results, and the need 
to restrict ourselves to values of In w < 3 — and a compromise was reached. This state of 
affairs was far from perfect, as the statistics generated from a mere 12 points were really 
not very convincing, and there were points towards the top left of the graphs for which 
lnu> was greater than 3 or so. Nevertheless, we present the results here for comparison 
with later results (Figs. [4.17] and 4.18|) . 

The graph for c appears to suggest c ~ 0(0.01) or thereabouts, which is much smaller 
than what we had before, from the MINIMAL results. The reason for this might be that 
12 points do not reliably fix the gradient on a graph. Note, also, that if we take c ~ 0.01, 
we have IncL rs 2.3, which requires Inw 2 -C 4.6. This was grossly violated in our graphs, 
as we allowed \nw 2 to go up to 6 or slightly beyond! Thus, these results are invalid and 
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Figure 4.17: Graph of c (derived from plots of \nw 2 vs. In J) against L 
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Figure 4.18: Graph of (derived from plots of lnu> 2 vs. In J) against L 
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Figure 4.19: Graph of Q against L, derived from plots of Jw 2 against Inw (the line at 
Q = 0.32 represents the theoretical value) 



useful information cannot be drawn from them. The statistics from these graphs actually 
point to their own invalidity. It would be impossible to restrict our results to Inw 2 <C 4.6 
as our results are all Inw 2 > 2 which already violates the regime. This forces us to take 
another approach. 

To overcome the restriction of the range of allowed Inw, we can plot Jw 2 against Inw 
for each value of L. This I have done, and the results are presented below: 



Gradient Intercept 



R 2 



1000 

900 

800 

700 

600 

500 



-0.16 ±0.02 
-0.15 ±0.04 
-0.15 ±0.04 
-0.14 ±0.10 
-0.18 ±0.10 
-0.12 ±0.06 



1.01 ±0.13 
0.98 ±0.13 
1.01 ±0.12 
1.08 ±0.28 
1.16 ±0.26 
0.92 ±0.14 



0.559 
0.348 
0.383 
0.0799 
0.112 
0.144 



^From Eq. (|4.7p , f2 should be 1/ir « 0.32. We plot the gradient as a function of L in 
Fig. |4.19[ There are three sets of points: omega, min omega and max omega. The min 
and max series represent the minimum and maximum values omega can take, within the 
bounds of experimental error. Thus, they can be thought of as the ends of (imaginary) 
error bars. The horizontal line in the middle of the graph at Q = —0.32 represents the 
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Figure 4.20: Graph of c against L, derived from plots of Jw 2 against Inw 



theoretical value. One can see that all the results agree with this predicted value, albeit 
with rather poor precision. Coupled with the earlier results from MINIMAL, this is mildly 
encouraging. 

^From the intercepts, we can work out c. We may, with some confidence, take Q to 



be 1/tt. The values of c obtained thus were plotted against L, in Fig. |4.20| . Because of 
the exponentiation of the intercept involved in obtaining c, errors in the intercept were 
greatly magnified, leading to relatively poor precision in our determination of c. However, 
it agrees with the earlier result of c ~ 0.5, though this set of results seems to point, rather 
vaguely, to c ~ 0(1). 



4.12 MINIMAL2: Cornell's approach 

S. Cornell, in his paper ||, performed extensive simulations of diffusion-annihilation re- 
actions, and described in detail his method of tackling the problem. His method, which 
allows particles of both species to diffuse over the same time steps, can lead to crossed- 
over buyers and sellers, which need then to be removed, along with the usual overlapping 
traders. One possible advantage of such a method, compared with MINIMAL and MIN- 
IMAL 1, is that the asymmetry between buyers and sellers, in which one type of trader 
always diffuses before the other type, is removed, and both are allowed to diffuse simul- 
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taneously. However, an inevitable consequence of this, when used with discrete time, is 
the crossing over of traders. At first, one might be tempted to think that such a scenario 
cannot occur in the market because traders generally try to sell as high as possible, and 
buyers buy as low as possible. However, because exchange of information via the trad- 
ing screen is not instantaneous (the screen is updated every so often), and humans can 
sometimes make mistakes, it is possible to have the best bid at a higher price than the 
best offer. In this case, usually the cross-over would be spotted by a broker, who would 
then inform the two traders of this situation so that they can do a deal immediately. The 
crossing-over characteristic of Cornell's model, which is perfectly legitimate in chemical 
reactions where annihilation is not immediate but follows an exponential decay law, is 
actually not as far-fetched as it might seem, when applied to financial markets. 

Strictly speaking, this is not the model described in flQ |, but a quick modification 
of the program was performed to see what such a method would yield. After several 
runs of this method the results did not look much different from those obtained using 
MINIMAL1. Further investigation is required before any conclusions may be drawn. 



4.13 The scaling law for NUM 



One interesting scaling law not mentioned in [K| concerns the equilibrium number of 
traders in the market, or NUM, where it is understood to mean the value at equilibrium. 



This is clearly shown in the two graphs, Figs. 4.21 and 4.22. Fig. 4.21 shows that 



NUM oc J and the gradient is the coefficient of proportionality. From plotting many 
such graphs of NUM against J for different L, one can collect the gradients (NUM/ J) 
together and thus plot ln(NUM/ J) against InL (Fig. f4.22|) . The graph is a straight line, 
with gradient 1.999±0.007 and intercept -0.688±0.004. We may infer that NUM/ J oc L 2 
with a constant of proportionality of 1/2 (from exp(— 0.688)). Note the remarkably good 
correlations on the graphs and the negligible experimental errors in the determination of 
the gradients. We may summarize, with great confidence, that 

JL 2 



NUM 



(4.{ 



obtained empirically. 

Dimensional analysis tells us that NUM is dimensionless. The right hand side of 
Eq. (|4.8|) , however, is not dimensionless. We can make it dimensionless by introducing D 
in the denominator, as D/J and L 2 have the same dimensions. The diffusion coefficient 
that we have been using is D — 1/2, so the true scaling law is 

JL 2 



NUM 



AD 



(4.9) 
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Figure 4.21: Graph of NUM against J for L = 1000 
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Figure 4.22: Graph of \n(NUM/ J) against InL 
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There is a firm theoretical justification for this scaling law. Consider a trader entering 
the system at the edge of price space. He will remain in price space until he hits the 
reaction front and annihilates with a trader of the opposite species. The reaction front 
is approximately in the middle of price space, a distance L/2 from the newly-inserted 
trader. From Einstein's work on Brownian motion (see jl3)), the mean square distance 
A 2 diffused by a particle in time t is ~ Dt. The time taken T for the new trader to diffuse 
to the reaction front is given by 



Now, we multiply both sides by J, the rate of insertion of traders, and divide by D: 



JT is simply the number of traders inserted between the time of the entrance of our trader 
into the market, and the time of his disappearance through annihilation. This is precisely 
the equilibrium number of traders in the market. One can see this in the following way: 
the total number of traders keeps increasing with every new trader inserted, up till a time 
T after the insertion of the first trader, when annihilation starts to occur, and balances 
the flux of new traders. This is only approximate, of course, to within a factor of order 
1. Thus, we have shown that NUM ~ JL 2 /AD. 

4.14 Scaling law for 75 

Now, we turn our attention to the scaling law obeyed by the instantaneous dealing time, 
ts, defined as the time between two consecutive annihilations. The program was modified 
so that it could record values of Ts and calculate their average and variance. There is a 
vital difference between this quantity and those we have been investigating thus far: the 
time between sales is not a function of the state of the system. All the other quantities, 
such as bid-offer spread and midmarket, could be read off the trading screen, and they 
were quantities that were characteristic of the state of the market. The instantaneous 
dealing time, however, is different, in that at any time one does not know when the next 
sale is going to be, and thus what the current value of t$ is. One can only record t$ as a 
historical quantity. Furthermore, the number of readings one can take of this quantity is 
not fixed, but varies, depending on J as well as the number of time steps the simulation 
is run for. 

The program was run for the range of values of J between 0.00004 and 0.007, which 
has been shown to produce w that satisfies 10 < w < 100, ensuring that the midmarket 
fluctuation, and other important lengths, are at least a factor of 10 or so away from the 
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size of the system, or the size of the smallest unit of length. For the smallest value of 
J, there were only 55 values of Ts obtained, whilst that number was 7140 for the largest 
value of J used. Therefore, the average and variance of t$ for the larger values of J were 
more reliable than those for smaller J. 

Originally, a graph of In Ts against In J was plotted for all the values, but while most 
of the points lay along a straight line, points corresponding to the 10 or so smallest 
values of J were rather scattered and did not lie so well on the line. This was because 
of the relatively small number of data points over which the averages and variances were 
obtained. These points were not as reliable as the others and so it was decided to exclude 
them in the graph plotting, and include only those for which J > 0.0002, all of which had 
been averaged over at least 210 points. A new graph was plotted. We obtained 



leading to the J exponent being —0.984 ± 0.005 and the constant in front of J being 
1.085 ± 0.035. The exponent is close to the expected value of -1, but unfortunately it is 
not quite within experimental error. The constant in front is also close to the expected 
value of 1, but just outside experimental error. We plot a graph of the fluctuation of time 
to midmarket sale, and find the following: 



The exponent is expected to be 2, since the fluctuation, defined as the variance of Ts, goes 
as the square of time, which is inversely proportional to J, from dimensional analysis. We 
obtain 0.71 ±0.08 for the constant in front, for which we have no "expected" value, since 
intuition does not tell us how scattered the times should be. Dimensional analysis is 
certainly correct in telling us the exponent of J, but it cannot give us any clue about 
what the constant should be, apart from a possible dimensionless parameter. 

It was noted that, because each of the simulations were run for a specified number of 
timesteps, namely 2 million, the total number of trades occurring within the time of each 
simulation was proportional to J. Thus, for smaller J, there were very few values of t$ for 
small J (as few as 55 for J = 0.00004), and more for larger J (7140 for J = 0.007). The 
number of readings from which each average and variance is obtained should not affect 
the accuracy of the results, though it might lower its precision, i.e. more scatter about 
the "true" mean and variance. In order to overcome the problem of having a variable 
number of Ts data points, the program was further modified so that the simulation would 
continue until a fixed number of trades had taken place. The limit was set at 1000, and 
the simulations terminated when 1000 trades had taken place and all their accompanying 
Ts values were recorded. The total running time was somewhere in the region of 18 hours. 



In t s = (0.08 ± 0.03) - (0.984 ± 0.005) In J 



(4.10) 



lnvar(r 5 ) = (-0.34 ±0.1) - (1.97 ± 0.02) In J. 



(4.11) 
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Figure 4.23: Graph of In r$ against In J for L = 1000 



The results obtained were plotted on two graphs, one of In ts against In J (Fig. 4.23| ) 



and the other of lnvar(rs') against In J (Fig. |4.24| ). When plotting the graphs, it was 
found that the data for the smaller values of J, though lying approximately along a 
straight line joining them and the other points, were rather more scattered about that 
line than the other points were. They appeared less reliable and were therefore excluded 
from the two graphs presented. Only values for which J > 0.0002 were plotted. From the 
gradients and intercepts, the following information was gleaned: 

lnr s = (0.08 ± 0.03) - (0.984 ± 0.004) In J (4.12) 
lnvar(r s ) = (-0.32 ± 0.07) - (1.97 ± 0.01) In J (4.13) 

In the case of Ts, the constant in front is 1.086 ± 0.03, whilst the corresponding constant 
is 0.73 ± 0.05 for var(rs). One can see a remarkable similarity between these results in 
Eqs. (|4.12| ) and (|4.13| ), and the earlier ones in Eqs. fl410D and (PH) . In fact, the only 



real difference is the reduction in the experimental uncertainty in the parameters. Thus, 
the second, more detailed, set of results not only confirms the validity of the first, but 
adds to it by narrowing the margins of uncertainty. 

It is fair, in the light of such compelling evidence, to believe that Ts goes as 1/J, as 
both the exponent and the constant were found to be very nearly -1 and 1, respectively. 
Such a conclusion is supported also by intuition, which tells us that rg = 1/J. The 
constant in front has also been confirmed. For var (ts), the exponent can be taken to be 
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Figure 4.24: Graph of lnvar(rs') against In J for L = 1000 

-2, while the constant in front is 0.73 ± 0.05. It seems very likely that neither quantity 
has a logarithmic correction. 



4.15 Scaling law for r rec i U ced 

There is another way of defining the instantaneous dealing time. It is identical to that 
given in §[4.14| for single annihilations, but differs for multi-annihilation. We can call this 



^reduced, a sort of "reduced" instantaneous dealing time. If the last annihilation contained 
n simultaneous trades, and the next one contained m simultaneous trades, then 

Treduced T ■ (^'l^) 

n + m — 1 

n and m may be considered the respective degeneracies of the trades. It would be inter- 
esting to see what statistics this new quantity obeys. 

We repeated the simulations for this reduced r, doing them for exactly the same range 
of J as we did for normal r. 1000 data points were obtained for each J, making each 
point as reliable as any other. A graph of lnr rc d U ccd against In J was plotted (Fig. |4.25| ). 



This graph is interesting in that the points together form a precise curve rather than a 
line. There is very little scatter, except in those points with the smallest J's. The rest 
of the points trace out a curve that tends to bend upwards as J is increased. If we only 
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Figure 4.25: Graph of In r re d U ced against In J for L = 



1000 



take the points with the smallest J in isolation, we obtain a gradient of -0.84. If we take 
the points with the highest J in isolation, we obtain -0.58 instead. Thus, there is some 
variation of the exponent with J. If by brute force we calculate a line of best fit, its 
statistics would be 

In r rcduccd = (-0.793 ± 0.009) In J + (1.71 ± 0.07) 

which gives an exponent of -0.793 and a "constant" of proportionality of 5.5 ± 0.4. 

Similarly, with the fluctuation (i.e. variance) in r red uced, the graph is not a straight line 



but a curve that curves upwards as J is increased. This is shown in Fig. [4.26| . Its gradient 
varies from -0.47 at its absolute minimum to -1.73 at its maximum. The mean gradient 
and intercept are 

lnvar(r rcduccd ) = (-1.32 ± 0.04) In J + (5.3 ± 0.3) 

giving a very imprecise estimate of the constant as 210 ± 70. The explanation of why the 
points do not fit a straight line well has not been found, and further investigations are 
needed before definite conclusions can be drawn. We will not pursue this question here. 
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Figure 4.26: Graph of In var(r rc d U ccd) against In J for L = 
4.16 The minimal model with asymmetric fluxes 



1000 



It was now decided to generalize the situation to one in which asymmetric fluxes of 
buyers and sellers are permitted. Intuitively, one would expect that, in the steady-state, 
the midmarket moves linearly with time towards one end or the other, away from the 
edge with the higher flux. This is a particularly productive phenomenon to investigate, 
since it will tie in very well with our later work on the two-liquid model, where different 
fluxes of limit-order and market-order traders will be introduced, in order to reproduce 
the well-known widening of the bid-offer spread in the prelude to a crash. The analytical 
result, Eq. ( |3.5| ), for the speed of the moving midmarket, is suggested in |L0| . This is the 
result we shall test. 

Initially, the program was run for very long times, in order to get a feel for what 
happened to the midmarket when the fluxes were asymmetric. It was observed that the 
midmarket did move, though it did not always do so at a constant speed. To illustrate this, 



two graphs of midmarket variation are shown: Fig. |4.27| and |4.28| . Fig. [4.27| appears to 
be quite straight and linear, but Fig. [4.28| does not. The latter exhibits some irregularity. 
In fact, there were many other examples of such graphs, where the midmarket movement 
was initially quite linear, but after some time, deviated from linearity. Most tended to 
bend downwards, indicating a higher speed of movement. As the scaling law Eq. (|3.5j) 
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gure 4.27: Midmarket variation with time for J = 0.001 and AJ = 0.0001 
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gure 4.28: Midmarket variation with time for J = 0.001 and AJ = 0.00001 
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Figure 4.29: Variation of best bid and best offer with time for J = 0.005 and A J = 0.0003 



applies only to small times, T < L 2 /D, we should restrict our attention to times less than 
2 million, and for these, the graphs were approximately linear. 

It is significant that the introduction of asymmetric fluxes of buyers and sellers has 
not caused any change in the bid-offer spread. A graph was plotted of the variation of the 
best bid and best offer with time (Fig. |4.29|) , from which it may be seen that the spread 
has not increased with time. In addition, a graph of spread itself was plotted as a function 
of time (Fig. |4.30j ). Spread appears to stay constant with time (after steady-state was 
reached) . 



4.16.1 AJ dependence 

Our first task was to investigate the dependence of the speed of midmarket movement on 
A J. This was done by fixing all the other variables (J and L) and allowing A J alone to 
vary. A graph would then be plotted of In \dx/dt\ against In A J, from which the exponent 
of the A J dependence could be deduced. A wide range of values for AJ were tried, but 
it was found that not all such values produced good results. For A J / J that was greater 
than 0.4 or so, it was seen that the midmarket moved so quickly that it hit rock bottom 
(price = 0) in very little time, sometimes less 1 million time steps, which would have been 
insufficient time for the system to come to any sort of steady-state. In the past, it was 
observed that 1 million time steps or so were required before the system could come to 
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Figure 4.30: Variation of bid-offer spread with time for J = 0.005 and A J = 0.0003 

equilibrium; now, although equilibrium can never be reached in this asymmetric situation, 
a sort of steady state can reasonably be achieved in that time. At the other end of the 
scale, it is possible to choose a AJ to J ratio that is too small. The midmarket drifts 
so slowly that the gradient of the graph is negligibly small. Although not a problem in 
itself, when one plots the logarithm of this quantity, it is a large negative number. A 
small uncertainty in the gradient is magnified into a large uncertainty in the logarithm. 
In practice, what this means is that the points for which the ratio is smaller than 0.02 or 
so are quite scattered and do not provide much useful information. Therefore, we shall 
limit our investigations to the range 0.02 < A J/ J < 0.4. 

Our choice of J is still subject to the restrictions described in the earlier sections, that 
the lengths produced by such a J (e.g. w, Spr) must be much larger than a, the lattice 
spacing, and much smaller than L. Thus, we require 0.007 < J < 0.00004. In order to 
allow AJ to range over as large an interval as possible, it is wise to choose a large J, as 
AJ is limited by it. It was decided to use J = 0.005, and L = 1000, as usual. AJ was 



allowed to range from 0.0001 to 0.002. The graph is shown in Fig. |4.31| . As one can see, 
there is noticeably more scatter in the points towards the lower end of the graph (small 
A J). The statistics from the graph indicate that 

In \dx/dt\ = (0.98 ± 0.05) In AJ - (1.7 ± 0.4). 



Thus, the observed exponent of A J is in accordance with the theoretical value of 1. The 
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Figure 4.31: Graph of In \dx/dt\ against In A J for J = 0.005 



constant of proportionality appears to be 0.18, with the minimum being 0.12 and the 
maximum 0.27. 

There is no reason why we should not try to get more accurate and precise statistics 
from our results, which are perfectly valid. Looking at the graph, one can see that apart 
from the first 10 points or so, the data do fit a line very well, and that the fit is spoilt 
by the scatter in the lower 10 points. Therefore, we plot another graph, Fig. |4.32| , of the 
same two quantities, this time without the first 10 points. The statistics are: 



In \dx/dt\ = (0.99 ± 0.05) In AJ - (1.61 ± 0.35). 

The exponent is now even closer to the expected value of 1. The constant in front is 0.20, 
with the limits being 0.14 and 0.29. In the light of such compelling evidence, we can 
conclude that there is no logarithmic dependence of any sort on A J and that \dx/dt\ is 
proportional to A J. Moreover, we can now plot a graph of \dx/dt\ against A J directly, 
from which we hope to obtain a more precise estimate of the constant of proportionality. 
This was done in Fig. |4.33[ Because we are no longer plotting the logarithm of small 
values of \dx/dt\, it is possible to re-include the first 10 points in our graph. The graph 
gives us an intercept of 2.4 x 10~ 6 with an uncertainty of 4.4 x 10~ 6 , which makes the 
intercept effectively zero. The gradient is 0.209 ± 0.005. Thus, plotting a direct graph 
can give us a much more precise estimate of the constant of proportionality. 
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Figure 4.32: Graph of In \dx/dt\ against In A J for J = 0.005, without the first 10 points 
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Figure 4.34: Graph of In \dx/dt\ against In J for A J = 0.0001 



4.16.2 J dependence 

Our next task is to determine the J dependence of the formula for the speed of a moving 
midmarket. An inverse dependence is expected. The limits on the ratio AJ to J dictate 
that the maximum value of J used be J ma x — 50AJ and the minimum be J m in — 2.5AJ. 
We require that J ma x < 0.007 in order to keep all lengths at least a factor of 10 away 
from L or a, implying A J < 0.00014. Thus, in the interests of maximizing the range of J 
available to us, whilst using a "round" number, it was decided to set AJ = 0.0001, with 
J min = 0.00025 and J mnm = 0.005. 



The results were plotted on a graph of \n\dx/dt\ against In J (Fig. |4.34| ). This gave 
us the following: 

In \dx/dt\ = (-1.01 ± 0.07) In J - (16.1 ± 0.5) 

from which a constant of proportionality of 9.9 x 10 -8 could be obtained (limits 5.9 x 10~ 8 
and 16.5 x 10~ 8 ). Once again, the exponent confirms the analytic result. In order to refine 
our estimate of the constant of proportionality (which will be of use later), we do a direct 
plot of \dx/dt\ against 1/J. Here, it was noticed that the 5 points corresponding to the 
smallest J values we used were very scattered on that plot, so it was decided to exclude 
them (Fig. |4.35| ). It does indeed yield a more precise estimate: (9.3 ± 0.9) x 10~ 8 . The 
intercept, found to be (0.5 ± 1.1) x 10~ 5 , is effectively zero. Therefore, we may conclude 
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Figure 4.35: Graph of \dx/dt\ against 1/J for AJ = 0.0001 (without first 5 points) 



with confidence that the speed is inversely proportional to J. 



4.16.3 L dependence 

Finally, we test the L dependence of the formula. This is, again, expected to be an 
inverse relationship. We tried many different values of L, from 100 up to 2000. Initially, 
J = 0.001 and A J = 0.00002 were used, but it was observed that the results were unusable 
for small L, such as 100 or 200, because there was far too much fluctuation from the mean 
midmarket. This is only to be expected from our earlier work on the fluctuation of the 
midmarket, w 2 , which is known to depend on 1/J. Therefore, it was decided to increase 
J whilst keeping the ratio A J / J the same in order to reduce the fluctuation from the 
mean midmarket. 

J = 0.005 and AJ = 0.0002 were used, which allowed the ratio AJ to J to stay 
within the allowed limits (see § 4.16. 1| ) for 100 < L < 2000. A log-log graph was plotted 



(Fig. [4 . 3 6|) . The line is fairly straight, but with some scatter, especially for small L. The 
reason is that for these small L, the midmarket moved towards zero so quickly that not 
much time elapsed before it crashed into the lower boundary. There were not many points 
from which one could derive a gradient, and in any case, there was much fluctuation 
in the instantaneous speed of the moving midmarket, and sometimes the fluctuations 
dominated, leaving the drift barely observable. All these factors made it very hard to 
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Figure 4.36: Graph of In \dx/dt\ against InL for A J = 0.005 



determine \dx/dt\. The gradient of the log-log plot is -1.2, which is not too bad, though 
we can certainly do better by excluding from our graph the first 6 points with the smallest 
L, which are the most scattered. This amounted to ignoring all the results for L < 300. 
This new graph is shown in Fig. [4.37|. The improved statistics are: 



In \dx/dt\ = (-1.02 ± 0.06) InL — (3.0 ± 0.4). 

The exponent is now much closer to the theoretical value of -1. The constant in front is 
0.05 ± 0.02. A more precise estimate may be obtained by a simple plot of \dx/dt\ against 
1/L, as we discovered earlier. Once again, the 6 points corresponding to L < 300 (now 
appearing last) will be excluded from the plot, which is shown in Fig. [4.38| . The intercept 



is (—0.9 ± 7) x 10~ 6 which is basically equivalent to zero. The slope is 0.045 ± 0.004, thus 
giving us an estimate of the constant of proportionality that is 5 times more precise than 
that provided by the log- log plot. We are therefore satisfied that the L dependence really 
is an inverse one. 



4.16.4 Constant of proportionality 

Having established the functional dependence of the speed of the moving midmarket, 
\dx/dt\, on J, A J, and L, we have demonstrated the truth of the formula suggested by 
Kogan and Eliezer, to within a constant of order 1. We shall now attempt to determine 
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what that constant is. It is not necessary to calculate numerically the D dependence 
(which we have never done for any quantity) because dimensional analysis alone tells us 
that if the speed, with units of dollars/sec, is proportional to A J/ JL, there must also be 
a factor of D there, in order to balance the units. 

For each of the dependences we investigated, a partial constant of proportionality, 
which included in it the other two variables, was obtained. ^From each of these we may 
infer an estimate of the constant in front of all three variables, i.e. the constant k, defined 

kDAJ 
x (t) = --^t. 

^From the first set of results, concerning the A J dependence, we deduce that k = 2.09 ± 
0.05, whilst for the second and third data sets, we obtain 1.86 ± 0.18 and 2.25 ± 0.2, 
respectively. We summarize this with a table: 

Set k Min Max 



1st (AJ) 2.09 ±0.05 2.04 2.14 
2nd (J) 1.86 ±0.18 1.68 2.04 
3rd (L) 2.25 ±0.20 2.05 2.45 

Although not tremendously precise, the data do tell us that the constant is close to 2, 
which was the factor suggested in Eq. (3.37) of |TIJ. Our most precise set of results was 
the 1st set, without a doubt, because all the points lay on a straight line, for both the 
direct and the log-log plots, and it is reasonable that we should pay more attention to 
its conclusions than to those of the subsequent sets, for which certain wildly inaccurate 
points had to be excluded. We can therefore conclude that our numerical simulations 
favour a constant of about 2.1, with an error of 0.1 or so. 



in 



This is supported by the theory. The formula in |10[ was obtained from Eq. (3.36) 
T0[ ] by setting ((x,t), the density difference between buyers and sellers, to zero for 



the midmarket, and then solving the resulting quadratic in x with the stochastic part 
(the Fourier series with noise-dependent coefficients) neglected. Expanding the square 
root in the quadratic formula gives Eq. |3.5| . Of course, for the two-liquid model, when 
the flux of market-order traders exceeds that of limit-order traders, the system decouples 
into two reaction fronts, between limit and market order traders, with asymmetric fluxes. 
Here, the bid-offer spread widens with a speed equal to twice that calculated here for the 
minimal model, because there are two fronts. We will investigate the two-liquid model in 
the next section. 
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5 Computer simulations of the two-liquid model 



Having spent much of our time and effort investigating the minimal model, it is now 
appropriate to turn our attention to the two-liquid model. It is more complicated than 
the minimal model to investigate, because there are four kinds of traders: limit-order 
buyer, limit-order seller, market-order buyer and market-order seller, and for this reason 
it is less analytically tractable. Many of Cardy et. al.'s methods which worked well with 
the minimal model cannot be applied here. Thus, there is good reason to resort to 
numerical simulations to solve the two-liquid model. 

The intention of the two-liquid model is to provide a more realistic description of a 
financial market in the prelude to a crash. In particular, what motivated the inclusion of 
the "second liquid" was the desire to reproduce the well-known phenomenon, observed by 
market practitioners, of the widening of the bid-offer spread. It has been demonstrated 
(Fig. |4.29| and |4.30|) that the asymmetric minimal model, though able to produce a steadily 



moving trade price, does not exhibit bid-offer spread widening. Market-order traders, 
believed to be the missing element, were added to reproduce this phenomenon. The rules 
of trading can be represented by the following equations, if we take B and S as the 
limit-order traders, and B' and S' as the market-order ones: 

B + S -> 
B + S' -> 
B' + S -»■ 

Note that B' + S' — > does not occur, since market-order traders do not put up their 
prices on the trading screen and so they cannot see one another. In numerical simulation, 
the problem arises of which ones to annihilate when one has more than two types of 
traders at a reaction front. For example, if there are B and B' at the same point in price 
space as S, do we allow the limit-order traders to trade first, or let the trades take place 
randomly? In the interests of computational efficiency, it was decided to allow limit-order 
traders to trade first, and then let the remainding LO buyers trade with MO sellers. The 
three types of annihilations occur in the order shown above. Although this might not be 
representative of a real market, it should suffice for a first attempt at simulating one. 

There are so many things in this model worthy of investigation that our research could 
do no more than scratch the surface, since our time was limited. Nevertheless, we did 
manage to get some interesting results, which are described in the following sections. 
We shall not dwell on the intricacies of the computer program used to perform these 
simulations, since the principles used here are no different from those used for the minimal 



56 



model, and those we have already discussed in ample detail in 



4.1 



5.1 Variation of spread with fraction of market-order traders 

It was decided to measure the bid-offer spread as a function of the fraction of market-order 
traders. Because we are using a model in which only the flux of traders is controllable, 
and the total number of traders of either type fluctuates, we can only measure the fraction 
of market-order traders with respect to the relative fluxes of the two types of traders. In 
other words, the fraction of market-order traders (MO) is here defined as the MO flux 
divided by the total flux of both types of traders. Using this definition, it is possible to 
change the fraction of MO traders by changing the MO flux, while keeping the total flux 
the same. It will be seen that this definition of MO fraction is supported by theory (which 
does not support a definition based on trader number). 

As time progresses, we know that in cases where MO flux exceeds LO (limit-order) 
flux, the bid-offer spread widens. The problem of what value to take as the spread for a 
particular flux configuration thereby arises. It was decided to define the spread here as the 
average spread between t = 300, 000 and t = 330, 000. This interval was chosen because 
it was believed to be sufficient time for the system to come to a sort of steady-state, yet 
not too long so that the two separating reaction fronts crash into the boundaries of price 
space. The mean flux was chosen to be 0.0005, a relatively high flux, because midmarket 
fluctuations, w 2 , which we wanted to minimize in order to have more accurate results, 
are known to go as 1/J. Thus, a 10% MO flux would correspond to J mo — 0.0001 and 
Jlo — 0.0009. Simulations were performed for many different fractions and the results 
were plotted on a graph (Fig. |5J]) . 

The critical point seems to be approximately 0.5, above which the spread begins to 
grow very rapidly [j. There is some fluctuation in the results, but they do fit a curve quite 
well, though they would fit a straight line probably just as well — we cannot tell at this 
stage. The value of the spread is undefined for 100% market-order traders. 

It is interesting to consider what would happen to the spread as the MO fraction 
tends to 1. We may apply Eq. |3.5| to the two-liquid model if we assume that it consists 
of two decoupled minimal systems, of MO against LO traders, with reaction fronts which 
move apart (this only works if J mo > Jlo)- We must also assume that the two fronts 
act independently of one another. This is only an approximation and does not apply at 
the beginning of the simulation, when the system is still trying to reach a steady-state. 

^Preliminary results that bifurcation point in the Two-Liquid Model occurs at / = 1/2 was obtained 
by Adrian McGowan and one of us (I.I.K.) in the spring 1999. 
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Figure 5.1: Variation of bid-offer spread with fraction of market-order traders 



Eq. |3.5| tells us that for a system with asymmetric fluxes J + AJ and J — A J, the speed of 
movement of the midmarket (which is effectively the same as the rate of movement of the 
reaction front, since the spread stays constant) is 2D A J/ JL. In our case, J mo = J + AJ 
and Jlo = J — AJ. Making the necessary substitutions, and remembering that the 
fraction / of MO traders is defined / = J mo/ (Jaw + Jlo), we obtain 

d(Spr) _ AD[2f - 1) 
dt ~ L 

which tells us how the spread increases with time. Integrating, we have 

40(2/ - l)t 



Spr(t) 



So 



(5.1) 



where So is the "natural" spread which exists in the absence of the drift, at / = 0.5. The 
formula predicts, therefore, that the graph of spread against fraction / would be linear 
for / > 0.5, with a gradient of 8Dt/L. In our case, this would have been 1260, if we use 
t = 315, 000. The spread at / = 0.5, So, is about 84, leading to a prediction of 714 for 
the y-intercept. 



We have already observed that the portion of the graph (Fig. |5.1|) for / > 0.5 is more 
like a curve than a straight line. However, if we use a straight line as a first approximation, 
and fit it to the points / > 0.5, we can see how the theoretical prediction matches up 



with the data. This was done in Fig. 5.2 which was a plot of all the points for which 
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Figure 5.2: Variation of bid-offer spread with fraction of market-order traders, for / > 0.5 
and excluding 3 outliers 

/ > 0.5, excluding the last three points which were obvious outliers. The graph gives 
a gradient of 1354 ± 64 and an intercept of 654 ± 48. Comparing the gradient with the 
theoretical value of 1260, one can see that the agreement is quite good. The y-intercept, 
too, is in surprisingly good agreement with the calculated value. Having said that, it is 
clear that Eq. (|5.1|) gives little more than an estimate of what the spread is. In any case, 
the equation only works in the region / > 0.5. 

A market practitioner would want a better fit to the data than a crude straight line. 
To this end, we have fitted the curve (less the 3 outliers) to a cubic, as shown in Fig. 
The fit is remarkably good. It cannot be determined with any certainty at this stage 
whether the portion of / > 0.5 is a curve, like a cubic, or a straight line. 

One wonders whether, in spite of the excellent cubic fit performed above, a single 
description of the variation of spread, or indeed any quantity, can be valid for all < / < 1. 
A cause for concern lies in the very sharp transition from a stable system / < 0.5 to an 
unstable one / > 0.5, with / = 0.5 being the critical point. We expect 0.5 to be the 
critical point, because an infinitesimal increase of / beyond this point would produce an 
imbalance of MO and LO fluxes, which would, in our model, inevitably lead to a steady 
widening of the spread. Yet is this the case? We see a sharp rise in spread at around 0.5, 
but does a phase transition really occur at exactly 0.5? The spread plotted was only the 
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Figure 5.3: Cubic fit to variation of spread with MO fraction, excluding 3 outliers 

average around t = 315, 000; perhaps we have been observing the spread too soon after 
the beginning of the simulation, and the system has not had time to reach a steady-state? 
These are all valid questions, which we now try to answer by performing more extensive 
simulations. 

5.2 More definite results for spread 

Admittedly, t = 300, 000 to t — 315, 000 is rather soon to start observing market be- 
haviour; we know that for symmetric fluxes, the market requires approximately 1 million 
time steps to come into equilibrium. We were forced to do so in order to maximize the 
range of /'s accessible to us. It was necessary to compare spreads measured at the same 
point in time for each simulation. Had we chosen a later time to observe, e.g. 1 million 
steps, we would not have been able to observe the variation of spread with / for / > 0.75 
or so, for which the spread at t — 10 6 would have exceeded 10000. If we look at / = 0.0, 
and compare the spread obtained by averaging from t = 300, 000 to t — 330, 000, and 
that obtained by averaging from t = 10 6 to t = 2 x 10 6 , we see that they are 45.5 and 

II One might think that the problem can be solved by simply increasing L. Not only does this increase 
the size of price space, but it also slows down the widening process (speed of widening ~ 1/L) and overall, 
the time taken for the spread to reach L goes as L 2 . However, the time taken for the system to reach 
equilibrium as a whole goes as L 2 , cancelling exactly with the increase in time gained from increasing L. 



60 



■o 

(S 

a. 

tfi 




I I I I I 

0.0000 .i GD G 0.2G00 0.3000 0.4-000 0.5000 0.6000 0.7000 

MO frsctjon 

Figure 5.4: Variation of spread with MO fraction, with spread averaged over 1 million 
steps 



41.7, respectively. An acceptable difference, if we bear in mind the inherent fluctuations. 
However, for / = 0.50, the spreads are 84.0 and 127.5, respectively. They differ by a 
factor of 1.5. Therefore, it is probable that the system had not reached a steady-state by 
the 300,000th time step. 

In order to probe this further, the spread calculations performed earlier were repeated, 
this time taking the spread as the average spread between 1 million and 2 million time 



steps (Fig. |5.4|) . This is the most comprehensive set of results for two-liquid spread so 
far. Looking at the graph, one is immediately struck by how well the points / < 0.5 fit a 
straight line. The reason for this is that the average over 1 million time steps is less prone 
to fluctuation than that over a mere 30,000 steps. We are limited in what region of / we 
can probe, since the spread shoots upwards very quickly once / exceeds 0.5, and reaches 
1000 around 0.64 or so. One can also see that the graph divides neatly into two straight 
sections: / < 0.5 and / > 0.5. This is encouraging, as it accords with what we predicted 
using the formula for asymmetric fluxes applied to the two-liquid model, Eq. (|5.1|) . Once 
again, we plot the two sections, / < 0.5 and / > 0.5 separately to determine the gradient 
and intercept (Fig. and |5.6|) . 

We consider first the section corresponding to / < 0.3, which is the straightest part 
of the / < 0.5 half. The full plot in Fig. [5^1 shows an almost perfect straight line. The 
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Figure 5.5: Variation of spread with / for / < 0.3, with spread averaged over 1 million 
steps 
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expanded version in Fig. |5.5| supports this, with the fitted line achieving a correlation 
coefficient in excess of 0.94. The gradient and intercept are, respectively, 47.6 ± 2.2 and 
40.4 ± 0.4. Initially, attempts were made to include points up to / = 0.4 but it was found 
that the points did not lie so well on a straight line. At the moment, we do not have 
a theory that can explain this straight portion. Note that the spread at 0.5 cannot be 
worked out by extrapolation from this formula, since the spread starts to deviate from 
this line at / ~ 0.4. 

Turning to the / > 0.525 portion, we note first that, contrary to previous results, the 
points quite clearly fit a line. The gradient and intercept are 5927 ±223 and —2928 ±130, 
respectively. Theory predicts that, for t ~ 1,500,000, the gradient d(Spr)/df would be 
8Dt/L = 6000. The intercept is So — 4Dt/L, the first term of which has to experimentally 
determined. We find that the equation describing / > 0.5 is 

Spr(f) = (5927 ± 223)/ - (2928 ± 130). 

The gradient agrees with theory, to well within experimental error. From Fig. |5.4j , we 
find, by visual inspection, that the spread at the point of symmetry, / = 0.5, is 140 ± 20. 
Incidentally, neither straight line section, when extrapolated to / = 0.5, gives the correct 
result, since 0.5 is in the curved section joining the two straight sections, and is larger 
than either straight line formula would predict. Now that we know So, we also know what 
theory predicts for the intercept of the / > 0.525 graph: it is —2860 ±20 (the uncertainty 
comes from our experimental determination of So). The actual intercept agrees well with 
our prediction. Thus, we are pleased to see good agreement with theory. 

In the light of these much more detailed results, especially Fig. [5T4"1 , it appears that the 
spread starts to rise rapidly just before f = 0.5, because we have already observed that 
the spread at that point of symmetry deviates from either linear regime. We recall that, 
in our method of simulation, whenever there was an abundance of traders that needed 
to be annihilated at a point in price space, we decided to always allow the LO traders 
to annihilate first. Such a technique would, in general, cause the best LO bid/offer to 
be further apart than it would be if we had not made that imposition and instead had 
allowed the LO and MO traders to trade with equal probability. This is because the 
spread is determined by LO traders alone, and if LO traders are more likely than MO 
traders to get annihilated, it is reasonable that the best bids and offers would 'snap back' 
more often, resulting in an increased average spread. Immediately, we can see that this 
preference for LO-LO trades only affects those trades where both LO buyers and LO 
sellers are present, i.e. when the spread is zero. Therefore, for / > 0.5, where the spread 
widens at a constant rate, and apart from at the beginning, the LO traders are completely 
separated, and we would expect the LO-LO preference to make no difference, since the LO 
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traders on either side are no longer in contact. For / just below 0.5, we expect no steady 
drift of the best bid or offer due to asymmetric fluxes, and so the situation is less clear. 
During all those times when the best LO bid and offer are not in contact, the priority 
given to LO-LO trades makes no difference. When they are coincident, however, more of 
the LO traders are going to be annihilated than they would otherwise. As a result, in 
such an annihilation where traders of all types are present, the number of LO traders will 
be more greatly diminished than that of MO traders. To consider what happens to the 
LO best bid, envisage a situation where, at the midmarket, the number of LO buyers is 
less than the total number of LO/MO sellers, which in turn is less than the total number 
of LO/MO buyers. Following our scheme, LO-LO trades occur first, followed by LO-MO 
trades. All LO buyers are eliminated. The only way for the best bid not to change would 
be for the number of LO buyers to exceed the total number of sellers. Had we used a 
random scheme of mutual annihilation, however, since buyers outnumber sellers, there 
would be no reason why some LO buyers should not be left. Similarly, if we perform this 
analysis for sellers, we would find a similar prejudice against the best offer staying where 
it is. In both cases, the number of LO traders has to outnumber the total number of the 
opposite type of trader in order for the best bid/offer to remain unchanged. For a random 
annihilation scheme, the only condition is that one type of trader outnumber the other, 
and then there is the possibility of either the best bid or best offer to remain unaltered. 
Therefore, it can be seen that increased spread is an artefact of the our simulation. 

By itself, the above analysis does not explain why the critical point might occur below 
0.5. To do that would require one to know quantitatively how the LO-LO preference 
affects the spread. Alternatively, to obviate this problem, it might be useful to change 
the simulation algorithm so that trades are indeed random, and that no preference is 
given to any type of trade. Such a modification would greatly increase processing time 
though it would probably resolve our uncertainty. 

It is interesting to note the numerical prefactor that is defined as the ratio of the 
spread at the point of symmetry (/ = 0.5) to that for the minimal model (/ = 0.0). From 
the results we already have, we can immediately give a value for it. The former is 140 ±20 
and the latter 40.4 ± 0.4, both coming from the graph of spread in the region / < 0.5. 
The ratio is therefore 3.47 ± 0.50. 

6 Conclusion 

The numerical analysis performed in this paper confirmed that w 2 , Spr 2 and vax(Spr) are 
all proportional to D/J, as predicted by dimensional analysis. The results also confirmed 
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the existence of logarithmic correction for w 2 , but not for Spr or var(Spr). This is not 
surprising, because the midmarket is ergodic and, given enough time, it will cover the 
entire width of price space. Increasing L increases the range in price space over which the 
midmarket can roam; thus, the midmarket variance diverges as L — > oo. The constant 
in front of the D/J in the w 2 scaling law was confirmed to be l/n. The logarithmic 
parameter for w 2 is c ~ 0.56, 0.32 < c < 0.97, for L = 1000. It appears that c ~ 0(1) 
for 500 < L < 1000. Spr and vai(Spr) do not seem to have logarithmic corrections. This 
is because Spr (and var (Spr), which is derived from Spr) is a fundamentally different 
quantity from the midmarket variance. Intuitively, changing L should not change Spr, 
since it is a property of the reaction front. The constants in front of D/J for Spr and 
var(Spr) are 3.615 ±0.065 and 0.95 ±0.02, respectively, whilst the approximate 'constant' 
in front of D/J for w 2 is 1.9 ± 0.2. We can therefore observe a hierarchy in the three 
major lengths: var(Spr) < w 2 < Spr according to 0.95 < 1.9 < 3.6. It was discovered 
that NUM scaled as JL 2 /AD, and a theoretical justification was given. This law was 
the most exact of all the scaling laws that the data supported; there was almost perfect 
correlation for all graphs of NUM against J and L. Time to midmarket sale (t$) was 
investigated, and found to be equal to 1/ J, as expected. Its fluctuation also scaled as 1 / J, 
but with 0.73 ± 0.05 as the constant of proportionality. Neither quantity showed signs 
of having logarithmic corrections. A further quantity, T re d U ced, a sort of scaled time, was 
investigated. Its scaling law seemed a little less straightforward, as the points on the log- 
log plot quite clearly followed a curve rather than a line. Further investigation is needed 
to determine its exact nature. The minimal model was extended by allowing asymmetric 
fluxes of buyers and sellers. Simulations verified the A J, J and L dependences, and the 
constant of proportionality, of the formula (Eq. ( j3.5|) in this paper) predicted by Kogan 
et. al. in [IIJ for the speed of a moving midmarket in the presence of asymmetric fluxes. 

Investigations into the two-liquid model produced widening of the bid-offer spread in 
the prelude to a crash, when the flux of market order traders exceeded the flux of limit 
order traders. Application of Eq. (|3.5f) to the spread yielded approximate, though uncer- 
tain, agreement. A cubic fit to the spread as a function of MO fraction was performed. 
Further, more detailed, results confirmed that the graph of Spr against / was made up of 
two straight sections, in agreement with theory. The predicted gradient and intercept of 
the second straight section were corroborated by experiment. The formula for asymmetric 
fluxes was found to describe the bifurcated regime well, but sheds no light on the equilib- 
rium regime, where / < 0.5. The critical point was found to be approximately / = 0.5, 
or perhaps just before. The ratio of spread for / = 0.0 and / = 0.5 was calculated to be 
3.47 ±0.50. 
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The analysis presented in this paper has covered five major scaling laws described in 
1(J. However, there are many more to be tested, and so there is much scope for further 



research. Other scaling laws that might be investigated include density near the best 
bid/offer as a function of deal rate, and higher correlation functions describing equilibrium. 
We have only scratched the surface of the two-liquid model; it would be interesting to see 
how the other scaling laws are modified in the prelude to a crash. Finally, JlO describes 
a third model, the bias model, which attempts to incorporate herd, or crowd, effects. It 
does so by modifying the diffusion operators to contain a drift, which depends on whether 
the last trade has moved upwards or downwards. This modification models momentum 
trading, where the traders have a tendency to go with the crowd. Such a model is designed 
to develop instabilities, and the drift parameter may be considered to be in a meta-stable 
basin, buffeted back and forth by random market movements (diffusion), until the drift 
is so large that it is knocked out of its basin, into the surrounding unstable region. Once 
there, it gathers momentum by positive feedback, leading to a crash. In future work we 
shall investigate the bias model. 
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A Appendix 

A.l The master equation and the diffusion coefficient 

From the master equation for the trader diffusion process, we may derive the diffusion 
equation and thus an expression for the diffusion coefficient D. 

Let P(x\t) be the probability of having a particle at point x. Let p be the probability 
for a particle to hop away from a given point, in each time quantum r. Each hop is to 
a point a distance a away from the original point of the particle. Let the probability 
of hopping to the left and to the right be equal, i.e. p/2. Therefore, the probability of 
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moving to the left is p/2, moving to the right p/2, and staying put (1 —p). We write the 
master equation by considering the probability for the particle to be at x at a time t + r: 

P(x\t + r) = ^P{x + a\t) + ^P(x-a\t) + (l-p)P(x\t) 
P{x\t + t) - P(x\t) = ^[P(x + a\t) + P(x - a\t) - 2P{x\t)\ 

We expand either side in a Taylor series about t and x. The odd derivatives cancel out on 
the right hand side, leaving just the even derivatives, which add. For small a and r, we 
may neglect terms of second order or higher in r on the LHS and fourth order or higher 
in a on the RHS. The partial derivatives are all evaluated at x and t. 



dP(x\t) p 
T ^^ + - = 2 



d 2 P{x\t) (a 2 \ ^ d 2 P{x\t) (a 2 ^ 



dx 2 \ 2 J dx 2 \ 2 J 

dP(x\t) (a 2 p\ d 2 P(x\t) 
di ~ V 27 J dx 2 

This is the diffusion equation for P(x\t), with diffusion coefficient D = a 2 p/(2r). In our 
simulations, p = 1 so D = a 2 / (2r), or D — 1/2, for a — r — 1. 



A.2 Proof of (X) = (X) and var(X) = (X 2 ) - (X) 2 

If we have a set of N results {X{\ and we choose to record only the averages over n steps, 
Xj, of the j th set of n values, such that 

X 3 = - J2 X i, ^<3<N/n, 

Tl 

i=n(j-l)+l 

then we may obtain the average ((...)) over the entire set of data by 

1 , 1 N/n 



This we may rewrite as 



I j = l 



1 7V/n _ 



which is simply the mean of the set of recorded averages (n is a constant). Similarly, we 
may obtain the variance of an entire set of data from knowing the sum of the squares, 
and the mean, of sets of n values: 

var(x) = —-\nr 
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) 



2 



N/n ~ \ 

(x 2 ) - (x) 2 - 



N/n ) 



Thus, the variance of an entire set of data is equal to the mean of the averaged squares 
minus the square of the mean of the recorded averages. Armed with these relations, 
we may reconstruct the mean and variance of the whole set of data from the averaged 
recorded data. 

References 

[1] M. Araujo, H. Larralde, S. Havlin, and H. E. Stanley, Scaling anomalies in reaction 
front dynamics of confined systems, Phys. Rev. Lett. 71 (1993), 3592. 

[2] , Reply to Comment by S. Cornell, Phys. Rev. Lett. 75 (1995), 2251. 

[3] P. Bak, M. Paczuski, and M. Shubik, Price variations in a stock market with many 
agents, Physica A (1997), no. 246, 430. 

[4] G. T. Barkema, M. J. Howard, and J. L. Cardy, Reaction- diffusion front for A+B — > 
in one dimension, Phys. Rev. E 53 (1996), R2017. 

[5] J. -P. Bouchaud and M. Potters, Theories des Risques Financiers, Eyrolles, Alea- 
Saclay, 1997 

[6] J. L. Cardy, Field Theory and Non- Equilibrium Statistical Mechanics, Lectures pre- 
sented as part of the Troisieme Cycle de la Suisse Romande, Spring 1999, unpub- 
lished, http:/ /www-thphys. physics. ox.ac.uk/users/JohnCardy/home. html 

[7] S. Cornell, Comment on "Scaling Anomalies in Reaction Front Dynamics of Confined 
Systems", Phys. Rev. Lett. 75 (1995), 2250. 

[8] , Refined simulations of the reaction front for diffusion-limited two-species 

annihilation in one dimension, Phys. Rev. E 51 (1995), 4055. 

[9] S. Cornell and M. Droz, Steady-state reaction- diffusion front scaling for mA + nB — > 
[inert], Phys. Rev. Lett. 70 (1993), 3824. 

[10] D. Eliezer and I. I. Kogan, Scaling laws for the market micro structure of the inter- 
dealer broker markets, 

Oxford preprint OUTP-98-64P, cond-matt/9808240; 



68 



Capital Markets Abstracts: Market Micro-structure (WPS) Vol.2 No.3 March 1, 1999; 
http:/ /papers, ssrn. com/paper, taf? abstract- id= 147135. 

[11] J. D. Farmer, Market Force, Ecology and Evolution, Adap-Org preprint 9812005 

[12] L. Gain and Z. Racz, Properties of the reaction front in an A + B — > C type reaction- 
diffusion process, Phys. Rev. A 38 (1988), 3151. 

[13] C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the nat- 
ural sciences, 2nd ed., pp. 1-11, Springer, 1997. 

[14] M. Howard and J. L. Cardy, Fluctuation effects and multiscaling of the reaction- 
diffusion front for A + B^O, J. Phys. A 28 (1995), 3599. 

[15] B. P. Lee and J. L. Cardy, Scaling of reaction zones in the A + B — > diffusion-limited 
reaction, Phys. Rev. E 50 (1996), R3287. 

[16] R. N. Mantegna and H. E. Stanley An Introduction to Econophysics: Correlations 
and Complexity in Finance , Cambridge University Press, 2000. 

[17] D. C. Mattis and M. L. Glasser, The uses of quantum field theory in diffusion-limited 
reactions, Rev. Mod. Phys. 70 (1998), no. 3, 979. 

[18] A. A. Ovchinnikov, S. F. Timashev, and A. A. Belyy, Kinetics of diffusion controlled 
chemical processes, pp. 1-9, Nova Science Publishers, New York, 1989. 

[19] F. Reif, Fundamentals of statistical and thermal physics, pp. 583-585, McGraw-Hill, 
1965. 

[20] J. J. Sakurai, Modern quantum mechanics, pp. 10-23, 89-94, Addison- Wesley, 1994. 



69 



